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numerical results is non-monotonic. To leading order the total energy emitted in the merger phase 
scales like rf and the spin of the final black hole scales like rj, where r\ = q/(l + q) 2 is the symmetric 
mass ratio. We study the multipolar distribution of the radiation, finding that odd-Z multipoles are 
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q increases. We introduce and compare three different definitions for the ringdown starting time. 
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resolution-dependent time variations in the fitted parameters of the final black hole. By cross- 
correlating information from different multipoles we show that ringdown fits can be used to obtain 
precise estimates of the mass and spin of the final black hole, which are in remarkable agreement 
with energy and angular momentum balance calculations. 
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I. INTRODUCTION 

More than thirty years after the first numerical simulations of binary black hole dynamics, the numerical relativity 
community is finally ready to compare binary black hole simulations with experimental data. Thanks to a series 
of recent breakthroughs, long term evolutions of inspiralling binary black holes that last for more than one orbit 
have been obtained with several independent codes, and accurate gravitational wave signals have been computed 

QiiiaiaaiEnni. 

The use of numerical waveforms as templates for gravitational wave detection requires large-scale parameter studies, 
and correspondingly large computational resources. The main current technical problems in the field are the efficiency 
of the numerical simulations and the development of a "data analysis pipeline" , connecting numerical simulations with 
analytical calculations of the early inspiral and late ringdown phases, and (eventually) with gravitational wave searches 
in actual detector data. To build a common language between the numerical relativity and data analysis communities 
we must develop a deeper understanding of the physical content of the simulations using analytical techniques, such as 
Post-Newtonian (PN) theory and black hole perturbation theory. A better analytical understanding of the simulations 
is important for many reasons: 

(1) To determine which regions of the parameter space (mass, spin magnitude and inclination, orbital separa- 
tion, eccentricity...) must be explored by numerical simulations, and which regions can be covered by (say) 
analytically-inspired interpolations of the numerical waveforms. This would obviously save a significant amount 
of computing time. 

(2) To develop optimal strategies for the construction of detection templates, using a combination of numerical and 
analytical techniques. 

(3) To understand details of the non-linear physics encoded in the strong-field merger gravitational waveforms, and 
extract as much science as possible from a detection. 

In this paper we focus on point (3), and we try to develop a general framework to quantitatively compare analytical 
calculations of the inspiral and ringdown waveforms with the "full" waveforms produced by numerical simulations, 
extending from late inspiral through merger and ringdown. 

Our work can be considered an extension of the recent analysis by Buonanno, Cook and Pretorius henceforth 
BCP). BCP studied simulations of non-spinning, equal mass black hole binaries starting out at three different initial 
separations. In this work we examine a larger set of simulations performed using the Bam code [U, Q and the moving 
puncture method. We consider seven different mass ratios (q = M2/M1 ~ 1 to q ~ 4 in steps of ~ 0.5) with 
initial coordinate separation D ~ 7 M, roughly corresponding to ~ 2 orbits before merger. For each mass ratio, the 
simulations were carried out at three different resolutions. To explore the effect of initial separation on the physical 
parameters of the remnant, we also consider two runs at separation D ~ 8M (for q — 2 and q = 3), and one run at 
separation D ~ 10 M (for q = 1). We typically use an extraction radius r cxt = 30M, with the exception of the q = 1 
run with D ~ 10 M , in which case we extract gravitational waves at r cxt = 30M, 40M and 50M. 

SectionllTlcontains details of our numerical setup. In Section al Al we study in some detail a well-known issue with the 
extraction of gravitational waveforms from numerical simulations: the problem of fixing integration constants when 
we integrate the Weyl scalar "J 4 twice in time to obtain the gravitational wave amplitude h. Fixing the integration 
constants to zero produces a systematic drift in h and in its first time derivative. This drift is sometimes referred 
to in the literature as a "memory effect", but this is somewhat misleading. The so-called "memory effect" is really 
due to numerical errors, wrong initial conditions and limitations of wave extraction techniques, and it should not be 
confused with the Christodoulou memory, which is a true (if typically small) physical effect due to the non-linearity of 
general relativity [Hj]. We find that the extraction radius is critical to reduce the amplitude drift, and that resolution 
only seems to affect the drift for low-amplitude components of the wave. 

In Section [TTT1 we study the inspiral-merger transition. We start by projecting the 2.5PN gravitational wave am- 
plitude for quasi-circular, non-spinning binaries [l3l . [Til [l5| onto spin-weighted spherical harmonics. In this way 
we obtain the spin-weighted spherical harmonic components of the Weyl scalar as PN series in the binary's orbital 
frequency: ipi, m = ipi ,m{MQ). We refer to this analytical expression of the gravitational wave amplitudes as the 
Post-Newtonian Quasi- Circular (PNQC) approximation (see Section [ill Al for details, and Appendix [A"l for a complete 
list of all the multipolar components). 

The PNQC approximation can be used in two ways. First, given the orbital frequency evolution f2(i), we can 
compute (an approximation to) the multipolar components ipi tm . Conversely, given the modulus of the wave amplitude 
\ipl,m\(t)> we can numerically invert the PN expansions to obtain a PNQC estimate of the orbital frequency: £1 ~ 
wpnqc ■ In Sections IIIIBI we compare wpnqc with two alternative estimates of the orbital frequency, first introduced 
in BCP: u>Dm (an estimate obtained from the gravitational wave frequency) and to c (computed from the punctures' 
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coordinate motion). Using these three different estimates of the orbital frequency, we study the convergence of the 
PNQC approximation. We find that, as in the point particle case [l6j], the convergence of the PN series is not 
monotonic. We also study the effect of resolution and wave extraction on the agreement between PNQC results and 
numerical results. We find that low resolution increases numerical noise in the frequencies and amplitudes at late 
times. A small extraction radius produces systematic errors at large separations, where gravitational wavelengths are 
longer, but it does not sensibly affect the ringdown phase. 

In Section [ill CI we study in detail the total radiated energy E to t and the final angular momentum jfi n as functions 
of the mass ratio, providing fitting formulas for each of these quantities. We also compare the energy and angular 
momentum fluxes with their PNQC estimates, and we study (both analytically and numerically) the multipolar 
distribution of the radiation. To leading order, we find that E tot ~ rf and jfi n ~ 77, where 77 = q/(l + q) 2 is the 
so-called symmetric mass ratio, and we provide fitting formulas for these quantities. As predicted by the PNQC 
approximation, odd-Z multipoles of the radiation are suppressed in the equal mass limit. As the mass ratio increases, 
higher multipoles (with / > 2) carry a larger fraction of the total energy: for q > 2, I = 3 typically carries ~ 10% of 
the total energy (see Table |l] below) . 

In Section IIVI we turn our attention to the merger-ringdown transition. During ringdown the waveform can be 
described as a superposition of complex exponentials, the quasinormal modes (QNMs). In [17j | we argued that Prony 
methods (which are well-known in signal processing) are in many ways "optimal" methods to extract QNM frequencies 
from a numerical signal. After explaining our choice of the fitting window, in Section IIVBI we use Prony methods 
and standard, non-linear least-squares fits to look at the time dependence of the final black hole's parameters. We 
find resolution-dependent deviations in these parameters from the values predicted by linear black hole perturbation 
theory. These effects may be due to non-linearities and/or to rotational mode coupling, but at present we cannot 
exclude the possibility that they are, more trivially, an artifact of finite numerical resolution. In Section IIV CI we 
show that, by cross-correlating information from different multipolar components of the ringdown waves, we can find 
an empirical "best guess" for the optimal time to estimate the final black hole's mass and angular momentum. We 
argue that, because of the no-hair theorem, this best guess corresponds to the last time when the angular momenta 
(or masses) obtained by fitting the dominant multipoles agree with each other. In support of this argument, we also 
show that estimates of the mass and spin of the final black hole based on QNM fits are in remarkable agreement with 
wave extraction methods. 

Black hole QNMs do not form a complete set, and for this reason it is not possible to define unambiguously the 
beginning of the ringdown phase. In Section IIVDI we consider three different definitions of the ringdown starting 
time, two of which have already appeared in the literature (but not in the context of binary black hole simulations). 
The first definition is based on looking for the time at which a QNM expansion provides the best fit to the actual 
numerical waveform, in the sense of a suitably defined norm [18j . Unfortunately, when applied to our numerical 
waveforms, this method is not particularly useful. The reason is that the norm is quite flat (and even worse, has some 
oscillations) over a wide range of starting times around the minimum. A second, more useful definition looks for the 
time maximizing the energy content of the QNM component of the waveform. For this reason, following Nollert , 
we call it the Energy Maximized Orthogonal Projection, or EMOP. We find that the "EMOP time" £emop and the 
maximum fraction of energy carried by ringdown (~ 42%) are remarkably independent of the mass ratio q. This is an 
indication that the ringdown waveform is in some sense "universal" : it does not depend too much on the details of 
the pre-merger phase. To our knowledge, the third definition of the ringdown starting time has not been introduced 
before. It uses a detection-based criterion, maximizing the "effective energy" deposited in a matched filter. 

In the conclusions we present a list of open problems and directions for future research. 

To improve readability, some lengthy equations and technical material are presented in the Appendices. 

Appendix IA1 lists the spin- weighted spherical harmonic components of the Weyl scalar, up to and including 2.5PN 
terms in a PN expansion of the waveforms. 

Appendix[B]provides fits for the energy, angular momentum and linear momentum radiated after the estimated time 
of formation of a common apparent horizon (CAH). Since the total energy radiated in a simulation depends on the 
initial separation of the binary, in this Appendix we also try to provide estimates for the energy, angular momentum 
and linear momentum radiated "after plunge". A problem here is that the Innermost Stable Circular Orbit (ISCO) 
is a controversial concept for comparable-mass binaries, and there is no unique way to define the beginning of the 
plunge phase. Given these intrinsic ambiguities, we estimate the starting time of the plunge, iisco: as the time when 
the orbital frequency SI becomes larger than the ISCO frequency computed in PN theory (at 2PN or 3PN order, to 
bracket uncertainties). We also present a comparison of our results with PN estimates of the post-plunge radiation 
by Blanchet et al. f20j . 

Computational resources and resolution limitations reduce the accuracy of numerical simulations for large mass 
ratio. Unfortunately, many astrophysical black hole binaries could have q = 10 or larger (see eg. |21| and references 
therein). It is important to determine the maximum value of q that should be simulated in numerical relativity, or 
equivalently, the smallest value of q for which black hole perturbation theory can be considered adequate for detection 
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and/or parameter estimation. Appendix [Cl collects some results from perturbation theory that may be useful in this 
context. We point out that, for large mass ratio, our numerical simulations seem to be in reasonable agreement with 
perturbative calculations of particles plunging with large angular momentum into a Schwarzschild black hole. 

Finally, in Appendix[D]we introduce quantitative measures of the polarization state of the waveform. We show that 
the polarization of the wave (as viewed from the normal to the orbital plane) is circular for both inspiral and ringdown, 
with the exception of the unphysical portions of the wave: the initial data burst and the final, noise-dominated part 
of the ringdown waveform. 

In all of this paper we adopt geometrical units (c = G = 1). Unless otherwise indicated, physical quantities 
are usually normalized to the total Arnowitt-Deser-Misner (ADM) mass of the system M. The ADM masses of all 
configurations presented in this study have been calculated using Ansorg's [25| spectral solver for binary black hole 
puncture data. Due to the spectral accuracy and the compactification of the coordinates, which facilitates evaluation 
of the ADM mass at infinite radius, the uncertainties in this quantity are negligible relative to those arising out of 
the numerical time evolution. 



TABLE I: Summary of the main results of this paper (see text). To convert from radiated momenta to kick velocities in km s _1 , 
the numbers in this Table must be multiplied by c/10 4 ~ 30 km s _1 . Further details (including estimates of the uncertainties) 
are given in the bulk of the paper. 





jfin 


jQNM 


Btot 


(Vol = 2, 3) 


Jtot 


^EMOP 


(Vol = 2, 3) 


'^EMOP 


it'Peviop 


E ISCO 


(Vol = 2, 3) 


J ISCO 


10 4 Pisco 


^filter 


q 


M 




M 


Mi 


M 


M 




m 


M 


1.0 


0.689 


0.684 


0.0372 


(96.3,0.4) 


0.246 


0.0185 


(97.3,0.7) 


0.0700 





0.0266 


(99.1,0.5) 


0.1231 





0.028 


1.5 


0.665 


0.664 


0.0340 


(94.8,2.0) 


0.229 


0.0174 


(91.9,2.1) 


0.0676 


2.32 


0.0245 


(97.1,2.4) 


0.1165 


2.05 


0.026 


2.0 


0.626 


0.626 


0.0286 


(91.8,4.6) 


0.196 


0.0142 


(91.5,5.1) 


0.0565 


8.77 


0.0208 


(93.5,5.4) 


0.0985 


3.17 


0.021 


2.5 


0.584 


0.581 


0.0238 


(89.2,6.8) 


0.167 


0.0119 


(92.4,7.6) 


0.0480 


9.42 


0.0172 


(91.1,8.0) 


0.0850 


3.68 


0.018 


3.0 


0.543 


0.544 


0.0200 


(86.8,8.7) 


0.143 


0.0103 


(85.4,9.3) 


0.0438 


8.99 


0.0146 


(88.1,10.1) 


0.0737 


3.92 


0.015 


3.5 


0.506 


0.509 


0.0170 


(84.6,10.1) 


0.124 


0.0089 


(84.2,9.6) 


0.0387 


8.29 


0.0123 


(86.2,11.8) 


0.0645 


3.88 


0.012 


4.0 


0.474 


0.478 


0.0145 


(83.2,11.3) 


0.108 


0.0078 


(82.1,10.4) 


0.0345 


7.49 


0.0106 


(83.4,13.0) 


0.0568 


3.82 


0.011 



For reference, we find it useful to summarize some of our main results in Table [T] There we list, for each mass ratio: 

(1) the dimensionless angular momentum of the final black hole Jfi n /M| n as estimated from wave extraction methods 
(jfin) and from QNM fits (j'qnm); 

(2) the total energy and angular momentum radiated in each simulation (E tot /M , J to t/M 2 ); 

(3) the energy, angular momentum and linear momentum radiated in ringdown, where the ringdown starting time 
is chosen according to the EMOP criterion (Bemop/^i Jemop/M 2 , Pemop/M)', 

(4) the energy, angular momentum and linear momentum radiated after plunge, where the plunge is defined by the 
location of the 3PN ISCO (E ISCO /M, Ji SCO /M 2 , ^isco/M); 

(5) the effective fraction of energy detected by a ringdown filter (Eaiter/M ). 

The Table also shows the fraction of energy being radiated in the two dominant multipoles (I — 2, 3). 
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II. NUMERICAL SETUP 

The sequence of numerical simulations of unequal mass black hole binaries studied in this work has been obtained 
with the Bam code @ using the moving puncture method [f| H[. Specifically, we study here a subset of the sequence 
used in Ref. Q to determine the maximum recoil resulting from the inspiral of non-spinning black hole binaries. The 
Bam code has been described extensively in Ref. Q and further details of the numerical simulations of the unequal 
mass binaries are given in Ref. [|[ . Here we summarize the model parameters relevant for our present study. 

A sequence of quasi-circular initial data of non-spinning black hole binaries is determined by the initial coordinate 
separation D, the mass ratio q of the black holes, and the initial momenta Pi of each black hole. Approximate values 
of Pi appropriate for quasi-circular orbits were calculated using the 3PN-accurate expression given in Section VII of 
Ref. ■ For most of the models we consider in this work, the initial coordinate separation is D ~ 7 M (denoted by 
"D7"). The mass ratio is varied from q ~ 1.0 to q ~ 4.0 in steps of approximately 0.5. In order to assess the impact of 
larger initial separations on our results, we also construct models with larger initial separation: D ~ 10 M for q ~ 1.0, 
and D ~ 8M for q ~ 2.0 and q ~ 3.0. We will denote these models by D10 and D8, respectively. The complete set 
of models is summarized in Table [TTJ 

TABLE II: Summary of the main physical parameters for the series of simulations studied in this work, q denotes the mass 
ratio, D is the initial coordinate separation, and J the total angular momentum. We also list the simulation time at which 
the orbital frequency equals the orbital frequency at the 3PN Innermost Stable Circular Orbit or ISCO, £i|coi an estimate of 
the time at which a CAH forms, £cah 5 the time at which the energy flux has a maximum, £fl ux ; and the time at which the 
modulus of the I = m — 2 mode has a peak, t pcs ^. All quantities are normalized to the ADM mass M. The final column lists 
the number N of orbits until the estimated time of formation of the CAH. 



q 


D/M 


J/M 2 


t'mco/M tcAu/M t au JM t 


peak /M 


N 


1.00 


7.046 


0.8845 


211.9 


215.0 


231.8 


234.0 


1.94 


1.49 


7.044 


0.8494 


213.5 


218.2 


234.3 


236.4 


1.96 


1.99 


7.040 


0.7870 


211.6 


217.5 


233.6 


235.4 


1.93 


2.48 


7.036 


0.7232 


213.2 


221.0 


236.2 


238.1 


1.96 


2.97 


7.034 


0.6649 


213.2 


223.2 


237.8 


239.8 


1.98 


3.46 


7.030 


0.6132 


215.0 


226.8 


240.6 


242.5 


2.02 


3.95 


7.028 


0.5679 


216.7 


230.5 


243.0 


244.8 


2.06 


1.00 


10.104 


0.9826 


972.7 


973.7 


992.3 


994.5 


5.93 


1.99 


8.086 


0.8214 


438.0 


443.9 


459.2 


461.6 


3.44 


2.97 


8.038 


0.6865 


392.2 


402.9 


414.5 


418.1 


3.23 



All models have been evolved in time using a resolution of Mi/22.4 near the punctures, where M\ is the puncture 
mass of the smaller hole. The models starting from an initial separation D = 7 M have also been evolved using 
resolutions of Mi/25.6 and Mi/28.8. In the remainder of this work we will refer to these resolutions as low (LR), 
medium (MR) and high resolution (HR) . Gravitational waves have been extracted in the form of the Newman-Penrose 
scalar ^4. Unless specified otherwise, we use an extraction radius r oxt = 30 M. We decompose the resulting ^4 into 
modes by projection onto spherical harmonics of spin-weight s = —2 (see Ref. Q for conventions) according to 

00 1 

Mr$ 4 = MrJ2 -MmW , 4) 1>i ,m • (2.1) 

2=2 m=— I 

In Fig. [T] we show examples of the resulting modes by plotting \Re(Mr ifri <m )\, the modulus of the real part of 
the waveforms. Except for the spurious initial wave burst (visible up to about t = 50 M) and for the final, noisy 
signal following the ringdown phase, the imaginary part of ipi, m is related to the real part by a phase shift of 7r/2 
(see Appendix [D] for a more detailed discussion of the polarization of the waveforms) . The figures demonstrate that 
the I = 2, \m\ = 2 modes dominate the gravitational wave emission in all simulations. Contributions due to higher 
order modes become increasingly significant, though, as the mass ratio is increased. In all of our models we find the 
strongest contributions of higher-^ modes to result from m = ±Z. An exception to this rule is the equal mass limit 
(q = 1), where odd- to modes (including the I = m = 3 component) are suppressed. For this reason, in the left panel 
of Fig. Q] we only show modes with 1 — 2, 4. 

In Fig. [5] we plot the modulus \Mripi^ m \ of the amplitude of some of the dominant multipoles. Two features of 
this plot are worth stressing: (i) the I — m = 4 mode does not have a single, clear damping time in the ringdown 
phase (this is particularly evident for q = 2.0); (ii) the amplitude modulation visible in the inspiral phase is induced 
by some eccentricity in the initial data. This eccentricity seems to decrease during the evolution, but estimates of the 
eccentricity damping are beyond the scope of this paper. 
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600 

t/M 




600 



FIG. 1: |Re(AfrVi,m)| for q = 1.0 (left) and q = 2.0 (right), 
strongly suppressed, and we do not show it. 



For the equal mass (q = 1.0) binary the I = m = 3 component is 





500 



FIG. 2: \Mrtj)i im \ for different mass ratios. Each plot shows only some of the dominant components: I = m = 2, 3, 4 and 
(1 = 2, m = 1). The initial burst of radiation is induced by the initial data, and the wiggles at late times are due to numerical 
noise. 



The late-time, exponentially decaying portion of the waveforms is the ringdown phase. As the wave amplitude 
decreases, numerical noise gradually starts dominating the signal. In order to exclude this noisy part from the fitting 
of damped sinusoids in the modelling of the ringdown part, discussed in Section HVl below, we introduce a cutoff time 
beyond which we no longer use the waveforms. The practical criterion to choose this late-time cutoff will be discussed 
in more detail in Section [TV] 

In order to assess the uncertainties arising from the discretization of the Einstein equations, we performed a 
convergence analysis for mass ratios q = 1, 2, 3 and 4. We find our results to converge at second order, with the 
exception of the equal mass (q = 1) case, where we have fourth order convergence. Second order convergence is 
demonstrated in Fig. [3J where we show the real part of the I = rn = 2 and I = m = 4 components of the scaled 
Newman- Penrose scalar Mr ^4. Here we choose the worst case (mass ratio q = 4) but the convergence is still quite 
good, especially considering that we do not apply any time shift to the waveforms. We similarly find second order 
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FIG. 3: Convergence plots for Mrip22 (left) and MrijJa (right). These plots show the differences between runs at different 
resolutions (as indicated in the inset), scaled to be consistent with second-order accuracy. They refer to run D7, mass ratio 
q = 4 and r oxt = 30M. 



convergence for the radiated, energy and momenta (see [3| for further convergence plots). We are therefore able to 
apply Richardson extrapolation and use the difference between the values thus obtained and the high resolution 
numerical results as estimates for the uncertainties associated with the finite differencing of the equations. 

For part of our analysis, we find it helpful to have estimates of the merger time of the black hole binary. The 
most reliable estimate would be the formation of a CAH. In order to reduce the computational cost, however, all 
simulations have been performed without using an apparent horizon finder, so that we need to rely on alternative 
estimates. In the case of equal masses, we follow Ref. [j| and use the lapse function a to estimate the black hole 
merger time as the time when the a = 0.3 regions around each hole merge. Unfortunately, this criterion does not 
generalize straightforwardly to unequal mass binaries. For these cases, we instead locate the time when the ratio of 
the radial and tangential speeds of the punctures is equal to 0.3, which corresponds roughly to the time when the 
black holes reach the "light ring" in the effective-one-body model (23|. This is discussed further in Section Hill below. 



A. Memory effects in subdominant multipoles 



The effect of a gravitational wave on a detector in the far field of the source is best described in terms of the 
transverse-tracefree part of the metric. The two polarization states, h + and h x , of the gravitational wave are related 
to the curvature, expressed in terms of the complex Newman-Penrose scalar VE^, by 

# 4 = h + - i'hx ■ (2.2) 

Here h+(t— z) = h^J = —hyy, h x (t — z) = fi^y = for a wave propagating in the ^-direction. Note that different 
conventions (typically for the Newman-Penrose scalar) are used in the literature, correspondingly leading to different 
relations with h+ and h x . Ref. j24[, for example, has a factor 2 in their Eq. (5.3). 

Given the Newman-Penrose scalar ^4 for a particular mode, we thus have to integrate twice in time to obtain h + 
and hx and, in consequence, fix two constants of integration, which correspond to the values and time derivatives of 
h + and h x at the initial time as functions on the (celestial) sphere. Integrations in time over ^4 are also required to 
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compute the radiated energy, linear and angular momentum from the radiation content: 



= lim 



= — lim 



r 

16tt 
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dVL 
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<9 / V±dt 




oo •/ — OO / 



(2.3) 
(2.4) 
(2.5) 

(2.6) 



The definitions above are based on time integrals which start in the infinite past (at time t = — oo), and thus capture 
the complete gravitational wave signal. Starting the time integrations at t = — oo corresponds to the limit of infinite 
extraction radius on the initial time slice — the slice would then extend all the way to spatial infinity, no part of the 
waveform would be lost, and it would take an infinite time for the waves to reach the extraction sphere. With our 
current setup of the numerical codes this situation cannot be handled, and we work with finite extraction radii. The 
constants of integration would then correspond to the signal that has been lost. In order to accurately compute from 
the Newman-Penrose scalar the radiated energy and momenta and the gravitational wave strain required by data 
analysts, it is thus necessary to understand the influence of these constants of integration, and ideally, how to choose 
them correctly 1 . 

Naively setting the constants to zero typically leads to a non-zero value and slope of h+ and h x after the passage 
of the wave. This effect will in general have contributions from the signal that has been lost due to a finite extraction 
radius, from numerical error, and from the inherent ambiguities of the extraction procedure at finite radius. Further- 
more, a time independent gravitational wave memory effect is also possible and has been described in the literature 
[H, 03] ( see a ^ so |26j)- Apart from the effect due to an improper setting of the constants of integration, the other 
effects will accumulate over time, which may allow for some discrimination. 

While the time independent phenomenon is physically expected (although it should be small), the time-dependent 
drift phenomenon appears to be counter-intuitive and we expect /i+ and h x to settle down into a stationary state at 
late times. All unphysical effects in the non-zero value and slope of h + and h x after the passage of the wave should 
converge away with resolution and increasing extraction radius. A rigorous convergence test could attempt to identify 
a remaining physical gravitational wave memory effect (indepent of time after the passage of the wave) . Consistency 
with the physical situation requires that the slope of h + and h x after the passage of the wave converge away with 
resolution and increasing extraction radius. 



TABLE III: Average slope of the I = 2, m 
regression. 



2 component of h+ for different runs and mass ratios, obtained from linear 



Koppitz et al. recently argued that the choice of integration constants may be important for an accurate calculation of the recoil 
velocity of the final black hole |25ll . 
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q run 


Resolution 




10 b a+ 




ioV 




1.0 DIO 
1.0 DIO 
1.0 DIO 


HR 
HR 
HR 


30 
10 
50 


7.11 
2.24 
0.91 


57.6 
57.3 
56.9 


3.67 
2.37 
1.42 


0.033 
0.038 
0.036 


1.0 DIO 
1.0 DIO 
1.0 DIO 


MR 
MR 
MR 


30 
40 
50 


7.12 
2.27 
0.94 


57.7 
58.1 
58.8 


3.67 
2.38 
1.42 


0.033 
0.038 
0.036 


1.0 DIO 
1.0 DIO 


LR 

T R 

LR 


30 
i n 

50 


7.09 

9 9fi 

0.94 


57.4 

^7 Q 

58.8 


3.68 

Z.oo 

1.42 


0.033 

U.Uoo 

0.036 


2.0 D8 
2.0 D8 
2.0 D8 


LR 
LR 
LR 


20 
25 
30 


1.35 
0.56 
0.27 


2.16 
2.19 
2.19 


5.57 
4.57 
3.20 


0.022 
0.029 
0.029 


3.0 D8 
3.0 D8 
3.0 D8 


LR 
LR 
LR 


20 
25 
30 


0.86 
0.36 
0.17 


1.38 
1.41 
1.38 


3.33 
2.83 
2.00 


0.013 
0.018 
0.018 
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FIG. 4: Effect of changing the extraction radius on the "memory effect" when we do not apply corrections to the integration 
constants. These plots refer to the dominant multipole (I = m = 2) of run D10 with q = 1.0. 



We study this effect in more detail by considering the D10 run with q — 1.0. In Fig. 0] we plot the resulting I = 2, 
to = 2 contribution for h+ obtained at different extraction radii. The figure demonstrates two important features 
of this memory effect. First, the linear growth starts right at the beginning of the simulation, indicating that the 
memory effect is indeed essentially due to a non- vanishing constant of integration, or that possibly it is accumulated 
already in the early stages of the wave pulse (including the artificial burst of radiation). Second, the slope decreases 
significantly if we use larger extraction radii. 

We next apply a least-squares fit of a linear function f{t) = ag + a\t to h+ and h x resulting from the simulations of 
models q — 1.0, D10; q — 2.0, D8; and q = 3.0, D8. The resulting slopes are labelled as and of respectively in Table 
IIIII The table demonstrates in the case of the model with q = 1.0, D10 that the coefficients are essentially independent 
of the grid resolution. Columns 5 to 7 of the table indicate, however, that their dependence on the extraction radius 
can be rather well approximated by power laws: ~ r cx> anc ^ a i ~ r ex- This discrepancy between the + and x 
polarization modes is rather surprising, because the circular polarization of the waves implies that they differ merely 
by a phase shift of tt/2. The key observation in this context is that this simple relation between h+ and h x applies 
to the inspiral waveform but not to the spurious initial wave burst. This is explicitly shown in Appendix |Dj We thus 
conclude that the slope is a consequence of the omitted early wave signal (the constant of integration) or the initial 
wave burst. We emphasize, however, that the decreasing impact of the initial data pulse at larger radii is not due to 
it being dissipated away by numerical viscosity, which would have manifested itself in a resolution dependence of the 
slope. 
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FIG. 5: Effect of changing the extraction radius on the "memory effect" when we do not apply corrections to the integration 
constants. These plots refer to a small amplitude mode (I = m = 4) of run D10 with q = 1.0. The left panel shows the effect 
of changing the extraction radius at fixed resolution. In the right panel, we change the resolution at fixed extraction radius. 



The picture is somewhat more complicated in the case of the I = 4, m — 4 mode plotted in Fig. [5l where we do not 
only see a significant dependence of the memory effect on resolution, but also a non-linear trend of h +1 suggesting 
contributions from the ambiguities of wave extraction or numerical error. Consistent with the issue being due to 
the ambiguities of the wave extraction algorithm, we observe a significant decrease in the memory effect at larger 
extraction radius. In contrast to the I — 2, m — 2 case, we did not find a simple systematic dependence of the 
coefficients on extraction radius and resolution. 

We conclude this Section with a discussion of alternative choices for the integration constants. For each polarization 
and each mode in the spectral decompositions we work with (to take care of the angular dependence of these constants 
of integration) we can fix the integration constants by demanding that the time derivative of h + / x vanishes at late 
times. This can be achieved by matching to a ringdown signal, or more heuristically by subtracting the time average 
or time dependent polynomials, such as those obtained by the fitting processes mentioned above. The second constant 
of integration, which may have a contribution from a physical memory effect, is typically very small, and can be set to 
zero for many practical purposes if the wave extraction radius is not too small. From our observations, we conclude 
that using a sufficiently large extraction radius is certainly a highly recommended way of reducing spurious memory 
effects. 
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III. THE INSPIRAL-MERGER TRANSITION 



The parameters chosen in Section [IT] for the initial data do not give perfect quasi-circular (non-eccentric) orbits. 
This problem has been discussed in various papers P, [T3, H3, HH . A simple way to visualize the residual eccentricity 
of the binary's orbit is to compare the punctures' motion with predictions for circular, Newtonian orbits. At leading 
order, the quadrupole formula predicts that the orbital radius should evolve according to 29]: 

5 r 6 

where rj — MiM 2 /M 2 = q/(l + q) 2 is the symmetric mass ratio. From the relation between the orbital radius r and 
the (Keplerian) orbital frequency SI, M — S7 2 r 3 , we get the ratio of radial and tangential velocities in unequal mass, 
circular orbit binaries: 

This formula is of course a rough approximation, being based on the quadrupole formula and assuming a Keplerian 
orbit 0. In Fig. © we show the ratio of radial and tangential velocities v r /v t obtained from the punctures' motion and 
from the Newtonian quadrupole prediction. Curves labeled "Newtonian" are obtained by replacing the punctures' 
orbital frequency SI = lo c (see Section [ill B II below for details of the definition) in Eq. (|3.2p . Curves obtained from the 
actual puncture orbital motion clearly oscillate around the Newtonian circular value, mainly because of the non-zero 
orbital eccentricity. A similar effect was observed in Fig. 6 of BCP. 
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FIG. 6: Ratio of radial and tangential velocities for D8 runs and for two values of the mass ratio (q = 2.0 and q = 3.0). 

At early times in the evolution, say (t — £ pca k) ^ — 100M, the ratio |t! r /wt| < 0.05, and the orbit is (to a reasonably 
good approximation) quasi-circular. At later times v r /vt grows, as the motion turns from inspiral to plunge. Given 
the computational cost of implementing an apparent horizon finder during the evolution, we use the following rough 
criterion to locate the formation of the apparent horizon. In the effective-one-body model [23], the ratio between 
radial and tangential velocities is v r /v t ~ 0.3 (so that the motion is strongly "plunging") at the light ring r = 3M. 
Since the light ring should be close to the location where a CAH forms, we simply define the time of formation of a 
CAH tcAH as the point where the ratio v r /vt, as computed from the punctures' orbital motion, becomes larger than 
0.3 (see also the related discussion around Table Hl|) . Fig. [6] shows that v r /vt rises very steeply in this region, so we 
expect the error introduced by our rough approximation to be at most of order a few M. 



A. The Post-Newtonian quasi-circular approximation for the inspiral phase 



In this work we will perform extensive comparisons of numerical waveforms with the PN approximation. For this 
purpose it is useful to decompose the Weyl scalar ^4 in spin-weighted spherical harmonic components according 
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to Eq. dHJ). The V;,m's can be obtained by taking two time derivatives of the PN gravitational waveforms h+ x 
according to Eq. (|2.2[) . and then computing 

Mripi im =Mr J sin 6d8d<j) -2YC.JO ,0)# 4 = Mr J sin 6 d6 d<t> ^Y^O , </>) (h+ - ih^ . (3.3) 

The azimuthal dependence of the PN waveforms has the functional form (/ ftdt - 2MQ\aQ/ft ) - cj) [l3|. Thus, the 
expansion of the waveform h = h + — ih x in spin-weighted spherical harmonics has a time dependence of the form 
exp[— im(J fldt — 2M£l In fl/fl Q )]. For consistency, we use the same convention 2 on spin-weighted spherical harmonics 
as in Refs. 0, For the dominant, I = m = 2 component of the waveform we get: 



Mr(h+-ih x ) 2 , 2 = 8W|?7(Mfi) 2 / 3 
-107 + 34?? 



21 



1 + 557/ ; 2 107 (M») 2 /3 + 27r( MQ) - 2173 + 7 f 5 3 1 V 204V (Mn)^ 

e -im(/ndt-2MninfJ/n ) (3 4) 



7T + roj^J (Mfi) 5/3 



Here fio is an arbitrary constant [131 ] and the orbital angular velocity SI is a time dependent quantity, the 3.5PN 
expansion of which can be found (for example) in [30| . 

Kidder et al. [151 ] recently corrected an inconsistency in the derivation of radiation reaction terms of Ref. [301 ] . They 
also argued that radiation reaction terms are in fact negligible in the 2.5PN waveform, since they can be absorbed 
into a 5PN contribution to the orbital phase evolution. The constant w in Eq. (|3.4p depends on whether radiation 
reaction terms are absorbed into a redefinition of the phase. If we include radiation reaction terms in the waveform by 
using Eq. (27) of [l5| (as recommended by Kidder et al.) then w = —24. If instead we neglect the radiation reaction 
contribution by using their Eq. (32), we find w = —8/7. In the following we present analytical results including all 
known contributions to the waveform (including the nj-dependent 2.5PN terms). Given the ambiguity due to the 
inclusion of radiation reaction terms, we decided not to include 2.5PN contributions in our comparisons with the 
I = m = 2 numerical waveforms. 

As stated earlier, to compute the projection of ^4 onto spin- weighted spherical harmonics we must take the second 
time derivative of expressions like Eq. (|3.4|) above. Noticing that the logarithmic term in the phase is of 4PN order 
[l3(], we will simply neglect it when taking the derivative 3 . One can then show that, up to 2.5PN order, 

i>i,7n = -m 2 n 2 (h + - ih x )i itn , (3.5) 

the only exception to this rule being the 2.5PN contribution to the amplitude of the I = m = 2 component. 

In our calculation of the amplitudes we discard terms of order C(AfSl) 14 / 3 (i.e. we compute all terms in a 2.5PN 
expansion of the gravitational wave amplitude, as given in [H|])- We only list the positive-m components of the 
dominant multipolcs, since ncgativc-m components are obtained by the symmetry property 

fl n ,- m = (-l)% m . (3.6) 

The small mass ratio limit of these results was obtained by Poisson and by Tagoshi and Sasaki [3lj ] . For comparable 
mass ratios, we find that the amplitudes of the dominant components are: 



2 This definition does not include the Condon-Shortley phase, an extra factor of (— l) m which is needed for agreement with the usual 
definition of scalar (s = 0) spherical harmonics. 

3 This is consistent with the PN order considered here. While it would be preferable to keep the logarithmic term, it introduces an extra 
unknown constant which we choose not to worry about in the present work. 
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(3.7b) 



(3.7c) 



(3.7d) 



where 77 is the symmetric mass ratio and we defined the phase (f> as 

i$ — p im(f Sldt-2MCtln.il/Cto) 



(3.8) 



The complete expressions of all multipolar components are listed, for reference, in Appendix [A] The leading order 
term in (|3.7ap is nothing but the quadrupole approximation: see eg. Eq. (24) of BCP. As predictable from symmetry 
arguments the odd-m multipoles, being proportional to SM/M, are suppressed in the equal mass case. To recover 
the spin-weighted expansion of the waveform h + — ih x one only has to divide these expressions by —m 2 il 2 . The 
only exception to this rule is the I = m = 2 component: the term proportional to 112/5 in (|3.7ap is the lowest-order 
correction due to the fact that the orbital angular velocity il is in fact a time dependent quantity. 

Terms with I > 2 and higher-order PN corrections provide a strong consistency check on both the PN expansion 
and the numerical results. First, they tell us if the PN expansion is a good approximation for higher multipolar 
components of the radiation (I > 2). Secondly, they can be used to check convergence of the PN expansion for 
any (l,m). If the series is convergent, for example, going beyond the so-called "restricted PN approximation" (i.e., 
including higher powers of (Mn) 1 / 3 in the expansion for 7^2,2) should yield better agreement with the amplitude 
predicted by numerical simulations. 

Notice also that multipoles which are formally of higher PN order are not necessarily subdominant. For instance, 
1P2 ,1 is of order (Mil) 3 . Based on power counting, this term should be comparable to 7/^33 and larger than 7^4,4. 
However the amplitude of these terms is proportional to (8/3)-^7r/5 ~ 2.11 for (I — 2, m — 1), 27 ^ynjl ~ 44.31 for 
(/ = 3, m = 3) and (1024/9)^/77/7 ~ 76.22 for (1 = 4, m = 4), respectively. At the maximum orbital frequency we 
are interested in (the ISCO frequency, which is of order Mil = Milisco — 0.1), the (2 , 1) amplitude is much smaller 
than the (4,4) amplitude: 7/^,1 — 0.06 (SM/M) 7^4,4. For this reason, in the following we will limit consideration 
to terms with I = m = 2, 3, 4. Figure [2] shows that the dominance of these terms is quantitatively confirmed by 
numerical simulations of the inspiral-merger transition. 



B. Estimates of the binary's orbital frequency from numerical simulations 

In the following, we estimate the orbital frequency il of a binary at any given time by three different methods, that 
we list below. 

(1) Orbital frequency from the gravitational wave frequency: il ~ tODm 

This estimate of il is based on the observation that the gravitational wave frequency in a mode characterized 
by azimuthal number m is cjgw — fnil. In practice, the calculation can be carried out in two equivalent ways: 

(i) Decompose each mode into a real amplitude and a real phase, ipi jTO = At >m e^q>(i<pi , m ). Then compute: 
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91 ,m 

dt 



(3.9) 



(ii) Alternatively, observe that if some frequency dominates the Fourier expansion of a signal, this frequency 
can be estimated by computing 



1 

Wfl m = Im 

m 



V'l 



(3.10) 



The latter method was used also in BCP, and it relies on the (implicit) assumption that the modulus of the 
complex mode amplitudes tpi iTn changes slowly compared with their phase. However, we verified that methods 
(i) and (ii) yield results which are basically indistinguishable from each other. In the following, when we refer 
to LO£> m we always compute Eq. (|3.10p by finite differencing. In any case, we verified for all modes that using 
Eq. (|3.9p would not produce appreciable differences. 



(2) Orbital frequency from the coordinate orbital motion of the punctures: ft ~ to c 

The idea here is to convert each puncture's motion in the (x , y) plane into polar coordinates (r ,</>), then compute 
There are two problems with this estimate of the orbital frequency. The first is that this definition obviously 
depends on the choice of coordinates, and we expect it to get worse as we get closer to merger. In our particular 
set of coordinates, the puncture motion agrees better with other estimates of the orbital frequency for large mass 
ratio, when the system becomes more similar to a test particle moving in Schwarzschild. The second problem 
is that, to compare the puncture coordinate frequencies against the other two, we need to take into account the 
finite time it takes for the waves to reach the extraction sphere. We simply estimate this propagation time to 
be At ~ f C xti but since the propagation speed may differ from unity and waves are not emitted from the origin, 
this introduces a (small) additional uncertainty in the comparison. 

(3) Orbital frequency from the Post-Newtonian Quasi- Circular approximation: Q ~ wpnqc 

BCP made the remarkable observation that, even very close to merger, the / = \m\ = 2 modes of the inspiral 
waveform can be well approximated by the standard quadrupole formula for a Newtonian binary in circular 
orbit. They computed the leading-order term in Eq. (|3.7a[) : 

\8/3-=F2ifa(t)-&0 



Mr i 2 ,±2 = 32J ^(MOf^e^ww-^o; ; ( 3 n ) 

where <f)(t) is the accumulated phase of the orbit with respect to some initial phase (/>o Ignoring logarithmic 
terms, 



f n(t')dt'. (3.12) 

Jo 



BCP pointed out that this simple Newtonian Quasi- Circular (NQC) approximation can be used in two ways. 
First, given the orbital frequency evolution ft(t) we can compute (an approximation to) the wave amplitude. 
Conversely, given the modulus of the wave amplitude, we can estimate the orbital frequency Q — wnqc by 
inverting the modulus of Eq. (|3.11[) , and check whether c^nqc agrees with the estimates loom (computed from 
the gravitational wave frequency) and u c (computed from the punctures' coordinate motion). 

We will show below that using Eqs. l|3.7a[) - (|3.7d[) and, more generally, the expressions listed in Appendix [Al 
their observation can be extended to all multipolar components of the radiation. Our approximation improves 
on the simple NQC estimate in that we include all PN terms in the expansion up to 2.5PN, but it still assumes 
that the orbits are quasi-circular. For this reason, we call it a Post-Newtonian Quasi- Circular estimate of the 
frequency. 

The three different estimates are shown in Fig. [7J to be compared with Fig. 7 in BCP. For this plot we choose 
runs D8 (since the inspiral part lasts longer) and we compare two values of the mass ratio: q = 2.0 and q = 3.0. 
After the final black hole's formation (roughly, for t > £cah) the system is no longer a binary, and therefore different 
estimates of the orbital frequency disagree with each other. The most physically meaningful quantity after merger 
is the gravitational wave frequency u>Dm- This frequency levels off to a constant, which is roughly proportional to 
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FIG. 7: Orbital frequencies from the puncture motion (Q = u c ), from the gravitational wave frequency (fi = (VDm) and from 
the PNQC approximation (fi = cjpnqc)- Here wpnqc is computed by inverting Eqs. (|3.7a|l . (|3.7b[l and (|3.7c|) , and keeping 
only the leading order in the PN expansions on the right hand side. Times are measured starting from t pca ^, the peak of the 
I — m = 2 mode amplitude for the given mass ratio (see Table HT|| . Horizontal lines mark 2PN and 3PN estimates of the ISCO 
frequency (as listed in Tables IXl and IXlj) ; the vertical dashed line marks (an estimate of) the CAH formation. The plots refer 
to runs D8 and two different mass ratios (q = 2.0 on the left, and q — 3.0 on the right). 



the fundamental I — m quasinormal frequency of the final black hole (see Section IfVI for a detailed analysis of the 
merger-ringdown transition) . 

The puncture coordinate frequency lo c is a reliable estimate at early times and large separations, whereas LUDm is 
initially noisy, being contaminated by spurious initial data radiation or noise from boundary reflections. However, 
the puncture coordinates provide a bad estimate of the orbital frequency already ~ 30M before the peak of the 
radiation. In this sense, when we are close to merger our coordinates are not as good as the generalized harmonic 
coordinates used in BCP. Coordinate choices having such a big impact, some care is required if we want to attach 
physical meaning to quantities like lo c . For example, BCP use the "decoupling point" where lu c separates from Wrj m 
to mark the beginning of the merger phase. With our particular choice of coordinates this decoupling point would 
occur much earlier, clearly invalidating the estimate. We argue that comparing uoum and (our best PN guess for) 
wpnqc should provide a coordinate- independent, more reliable estimate of the decoupling point. 

In Fig. [7| the PNQC frequency wpnqc is computed by inverting Eqs. (|3.7p and keeping only the leading order. This 
simple leading-order approximation is in excellent agreement with the other estimates (lu c and uJDm) until ~ 20 M 
before the radiation peak. At this point the orbit transitions from inspiral to plunge, and we cannot expect the PN 
inspiral calculation to provide the correct orbital frequency anymore. 

We expect the transition from inspiral to plunge to happen, roughly, when the binary's orbital frequency crosses 
the ISCO frequency. To estimate the ISCO we look for extrema of the 2PN and 3PN Taylor expansions of the binding 
energy (see Appendix [Bl for numerical values of MfiiscO) an d Section III A of Ref. [32| for a discussion of this and 
alternative methods of estimating the ISCO). The corresponding estimates are marked by horizontal lines in Fig. [7] 
Around the transition region, wpnqc (which is computed assuming that the motion is a slow, quasi-circular inspiral) 
should deviate more and more from oj c and Wjj m . 

This statement is made more quantitative in Fig. [5J where we zoom in around the ISCO region. Thick lines (the 
actual gravitational wave frequencies of the system in a multipole with I = m, divided by m) and thin lines (the 
PNQC estimates) are almost parallel to each other before the ISCO, and they deviate significantly as the orbit crosses 
the ISCO. The pre-ISCO agreement is better, and the post-ISCO deviation larger, for the q = 3.0 binary than for the 
q = 2.0 binary. This is in agreement with the physical expectation that, as the binary's masses become comparable, 
the very notion of an ISCO becomes less and less significant: roughly speaking, the system cannot be described 
anymore by the simple-minded picture of a "small" particle orbiting a larger black hole. 

An interesting question is if the agreement between the PNQC frequency wpnqc and other estimates of the orbital 
frequency, namely lo c and uj) ra , improves if we include additional terms in Eqs. (|3.7a[) . (|3.7b|) and (|3.7c[) . In other 
words, can we use different estimates of the PNQC orbital frequency i^pnqc to estimate the convergence rate of the 
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FIG. 8: Orbital frequencies from the puncture motion (Q = u c ), from the gravitational wave frequency (SI = (VDm) and from 
the PNQC approximation (fi = <^pnqc) in the region around the ISCO. Here ijpnqc is computed inverting Eqs. (|3,7a|l . ()3.7b|) 
and (|3.7c|) . and keeping only the leading order in the PN expansions on the right hand side. 



PN approximation? Conversely, if we substitute u>Dm into Eqs. (|3.7p to compute some PN approximation to the 
amplitude, does the agreement with the modulus of the numerical amplitude \Mripi^ m \ get better as we increase the 
PN order? We address these issues below. 



1. Convergence of the Post-Newtonian quasi-circular approximation 

In Fig. [S] we show a simple "visual" convergence test of the PN approximation. We substitute u>Dm into Eqs. (|3.7a|) 
and (|3.7b[) and we compute successive PN approximations to the wave amplitude as functions of time; then we compare 
the results with the modulus of the actual numerical amplitude \Mrif>i m \. 

At early times we clearly see oscillations in the PN estimates of the amplitude, that damp away as the binary 
evolves. These oscillations are due to ujd2 being very noisy near the beginning of the simulation (compare the early 
portion of Fig. [7]), and they would not be present if we used as a reference oj c , which is much smoother at early times 4 . 
The first PN correction is seen to deviate significantly from all other PN approximations. This is a general feature of 
PN expansions. The poor convergence properties of the PN approximation have long been known in the point-particle 
limit [16j |. where exact results can be obtained by simply integrating the Zerilli equation. Fortunately, in our case 
higher-order PN expansions (of order higher than 1PN for I — m = 2) are reasonably consistent with each other. 

A comparison of the PNQC orbital frequencies computed at different PN orders is also instructive. Let us assume 
that, within the accuracy of our numerical simulations, Loom is a good representation of the "true" orbital frequency 
of the binary 5 . If by increasing the PN order we find that wpnqc gets closer and closer to oJDm, this would provide 
an indication that the PN expansion is converging to the actual solution of the full, nonlinear problem. 

Fig. 1101 shows the relative deviation between wpnqc> as computed by keeping more and more terms in Eqs. (|3.7p . 
and the supposedly more accurate orbital frequency u>D2- Once again, at early-times we see oscillations in the relative 
deviation, that damp away as the binary evolves. The magnitude of the relative deviation |(wpnqc — ^D2)I^>D2\ can 
be taken as an indicator of the accuracy of the PN approximation. These plots confirm, from a slightly different 



4 The reason why we do not choose u) c as a baseline for comparison is that this frequency significantly deviates from the others close to 
merger, where the comparison between PN theory and numerical simulations is most interesting. 

5 In practice, of course, this is only true in an approximate sense. For example, in the early inspiral u>u m is heavily contaminated by 
the initial data burst and boundary reflection noise; finite differencing effects will introduce errors in the calculation of ujDm given the 
computed waveform; and finally, the waveform itself is obtained at finite extraction radius, introducing additional uncertainties. What 
we are really assuming is that, taken together, all of these effects are smaller than the errors introduced by the PN approximation to 
general relativity. 
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FIG. 9: Amplitudes obtained by substituting uiDm into the PNQC equations are compared with the numerical amplitude. All 
plots refer to runs D8. Figures on the left refer to q = 2.0, those on the right to q — 3.0. The top row shows amplitudes for 
I = m = 2, the bottom row for I — m = 3. Vertical lines mark 2PN and 3PN estimates of the ISCO. 



perspective, the non-monotonic convergence of the PN series. After the transition from inspiral to plunge (very 
roughly corresponding to the vertical lines, marking the estimated location of the ISCO at 2PN and 3PN) the PNQC 
frequency, which is only valid for the inspiral phase, clearly decouples from ujd2, and the relative error becomes much 
larger. Perhaps in the future, as the accuracy of numerical simulations increases and higher-order PN calculations 
become available, it will be possible to use the change in slope of |(o>pnqc — ^D2)/^>D2\ to monitor the occurrence of 
an orbital instability (the "plunge phase" ) in full general relativity. 

We already pointed out that our assumed "exact" orbital frequency, uju2, is in practice affected by various sources 
of numerical error: finite-differencing errors, the finite extraction radius and the initial data burst all introduce 
uncertainties. To bracket these uncertainties, in Fig. [11] we choose one of our longest runs (D10 for q = 1.0) and we 
study the effect of resolution and extraction radius on (o>pnqc — ^D2)I^bi- 

Two remarkable features emerge from this plot. First of all, at lower resolution the "wiggles" induced by initial data 
are still visible at later times. The second effect is perhaps the most important for matching numerical waveforms to 
the PN approximation, and for building template banks for gravitational wave detection. We see that small extraction 
radii produce a systematic bias (i.e., a larger deviation o/wpnqc from ujr> m ) at large orbital separations. This effect 
is easily understandable: the typical wavelength in the "early" inspiral part is of order A ~ lOOilf , which is actually 
larger than the typical extraction radii used in the present simulations. Such small extraction radii inevitably produce 
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FIG. 10: Relative deviation between the orbital frequency obtained from the PNQC inspiral formulae and the "true" frequency 
of the signal ujdi- Plots refer to runs D8 with q — 2.0 (left), q = 3.0 (right). 









FIG. 11: The effect of changing extraction radius (top) and resolution (bottom) on the errors. These plots refer to mass ratio 
q — 1.0, run D10. The extraction radii (r cxt = 30, 40 and 50) and resolutions (HR, MR, LR stands for high, medium and low 
resolution, respectively) are indicated in the inset. 



a bias in the waveform. We observed a similar, and probably related problem in the context of what we called the 
memory effect (Section III A[) . 
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C. Radiated energy and angular momentum 



1. Total radiated energy 



The total radiated energy computed from wave-extraction methods, and the energy radiated into each multipole I, 
are listed in Table IIVI and plotted in Fig. [T3J The error estimates listed in the table are obtained from Richardson 
extrapolation as described in Sec. |TT] (cf. also Table [V] below, where we also list the radiated energies and the final 
angular momenta computed using QNM fits). 



TABLE IV: Total energy radiated in merger simulations of unequal mass black holes, and percentage of energy in the different 
multipoles (normalized to the total energy radiated in I = 2, . . . , 5). The numbers refer to high-resolution D7 runs, and error 
estimates are obtained using Richardson extrapolation. In parentheses we list numbers obtained eliminating the initial data 
burst (in practice, we remove all data for t < to = 75M). 
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FIG. 12: Left: total energy E to t/M radiated in the merger (including the initial data burst) fitted by Eq. (|3.13p and Eq. (|3.14[) . 
respectively. Right: numerical data for the energy Ei/M emitted into each multipole I are compared with the fitting functions 
(13.151) . Error bars are estimated by Richardson extrapolation. Notice the suppression of odd multipoles as q — > 1. 



Table IIVI and Fig. clearly illustrate that the relative contribution of higher multipoles becomes more relevant 
as the mass ratio increases. As expected from symmetry considerations (and from the calculations in Appendix 
|A"|) . odd values of m are suppressed as the mass ratio q — > 1. The dominant components for all mass ratios are 
(I, m) — (2, 2), (3, 3), (4, 4). We often observe a non-negligible contribution (partly due to spurious initial data 
radiation) also in (I, m) = (2, 1), (2, 0), (4, 0), (5, 5). The initial data radiation burst can be eliminated by starting 
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the integration of the energy flux after the initial burst has passed. In Table [TV] we decided, somewhat arbitrarily, to 
start the integration at t — 75M. Changes in the starting time have a marginal impact on the results: at the level of 
0.1% for the (2 , 2) modes, and of about a few percent for the weakest modes (which have higher errors anyway). 

In practical applications it may be useful to have some fitting formulas for the total energy radiated and for 
the contribution of different multipolcs. Since the energy is proportional to |^ , 4| 2 , and the I = m = 2 component 
ip2 ,2 ~ Tj dominates the radiation, we expect the total radiated energy to be roughly proportional to rj 2 (recall that 
the symmetric mass ratio r) = q/(l + q) 2 tends to 1/4 in the equal mass limit). Indeed, it turns out that the total 
radiated energy in the merger E to t is fitted extremely well (deviations from the data being < 4%) by the function 
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marginally improves the quality of the fit, bringing the agreement with the data to the level of ~ 1% (see the left 
panel of Fig. [T2"|) . 

The different multipolar components are slightly harder to fit. Since again we expect the energy in each component 
to be proportional to the square of the amplitudes, the even components with I = m should be proportional to 
[4g/(l + g) 2 ] , and the odd components should scale with [q(q — l)/(l + g) 3 ] . After some experimentation with 
including higher-order corrections in 77, we found that the following functions provide a satisfactory fit of the data: 
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For I = 2, we find the best-fit coefficients to be (ci = 0.024231 , c 2 = 0.012163); for I = 4, they are (ci = 
0.0010294, c 2 = -0.0005328). For I = 3 we find [d x = 0.00017, d 2 = 0.10509, d 3 = 0.13990), and for I = 5 
(di = 0. 000012, d 2 = 0. 011020, d 3 = 0.000463). Numerical data for different multipolar contributions and the 
corresponding fitting functions are shown in the right panel of Fig. 1121 



2. Final angular momentum 
A good fit to the final angular momentum, for small mass ratios, was found in Ref. Q: 

j~ 0.089 + 2.4 j . (3.16) 
(1 + 9) 

In this paper we compute the final angular momentum in two ways. One estimate, that we denote by jfi n , is obtained 
by subtracting the total radiated angular momentum (as computed by wave extraction) from the total angular 
momentum at the beginning of the simulation. A second estimate j'qnm is based on QNM fits, and will be described 
in detail in Section HVl below. The actual data, together with error estimates based on Richardson extrapolation, can 
be found in Table El We find our estimates to be accurate within a few percent for q > 3, and even more accurate as 
q—*l. We found that very good fits (accurate to within ~ 1% of the numerical data) are given by 

2 

JQNM ^ 3 - 352 7T^j2 - 2 - 461 JYTW ' (3 ' 17£l) 

2 

]fin ~ 3.272—^- 2.075-^. (3.17b) 
(1 + 9) (1 + 9) 



The quality of these fits, and the very good agreement between the two different estimates of the final angular 
momentum, is illustrated in Fig. 1131 
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FIG. 13: Angular momentum estimated from a QNM fit and from wave extraction, and corresponding fits. Error bars are 
estimated by Richardson extrapolation. 



The functional form of Eq. (]3 . 1 T[) can be justified by a simple physical argument. Consider an extreme-mass ratio 
binary with the smaller body orbiting near the ISCO. The orbital angular momentum at the ISCO, in the small mass 
ratio limit and for non-spinning bodies, is given by 
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(3.18) 



In Appendix [B] we show that the numerical results for the angular momentum radiated after the ISCO are well fitted 
by Eq. (ED , that we reproduce here: 
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Therefore the angular momentum of the final hole should be well described by 
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which is remarkably close to the best fits ([3, 17ft . 
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3. Energy and angular momentum fluxes 



The purpose of this Section is to compare analytical PN estimates of the energy and angular momentum flux for a 
quasi-circular, unequal mass inspiral against the corresponding numerical calculations. To begin with, we summarize 
how to compute these quantities in PN theory and in numerical relativity. 

The gravitational wave energy flux emitted by a binary moving along an adiabatic sequence of circular orbits is 
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currently known analytically through 3.5PN order for non-spinning BHs in circular orbits [30]. It reads 
^tPnqc = ^ 2(MO)1 o/3| 1+ ^_^_35^ (MO)2 /3 +47r(MO) 

/ 44711 9271 65 /3 + _ 583 \ ff(mn)B/8 
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^ + ^^ 2 )-(W /3 }, (3.21) 



where 75 is Euler's number and rj denotes, as usual, the symmetric mass ratio. The numerical energy flux can be 
obtained from the mode amplitudes Mr ?pi jTn (t) as 



dE 



= E f e, ; = E E ^ m = -T-Ei^.™wi 2 < ( 3 - 22 ) 



at lbn 

1=2 1=2 m= — l l,m 



where Di iTn (t) is a dimensionless first time integral of ipi m (t) defined by 



1 r* 

D l>m {t) = — dt'Mr^ l>m {t'). (3.23) 



MJ 

Each term in the sum (13.22j) represents the multipolar contribution of a different mode. A PN estimate for the 
flux in each (I , m) mode can be obtained by using the expansion of ^4 in spin-weighted spherical harmonics. Using 
the same approximation discussed in Section [ill Al namely neglecting the logarithmic term in the phase, we get the 
folllowing PN estimate for the (/ , m) component of the energy flux: 

F S° = ^ = I^l^f, (3.24) 

where we use the best available PN expansions of Mr ipi jm , as listed in AppendixjA] The 2.5PN term in the I = m = 2 
waveform suffers from the usual ambiguity in w related with the inclusion of radiation reaction terms [15J . For this 
reason we only consider I = m = 2 corrections up to and including 2PN terms. 
For quasi-circular orbits, the PN angular momentum flux is simply 

F PNQC = l^NQC (3 25) 

The numerical angular momentum flux Fj can be computed using Eq. (|2.5[) . 

In Fig. [T?] we compare the total energy flux (as computed from our numerical simulations) with the PNQC energy 
flux (|3.2ip evaluated at 2PN and 3.5PN. The 2PN, 3PN and 3.5PN expansions are very close to each other, and to 
improve readability we decided not to display the 3PN results 6 . At each different PN order we evaluate Eq. (|3 . 2 1 [) by 
using two different estimates of the orbital frequency (uJDm and lo c ). 

Some features of Fig. [TJ] should be quite familiar from the discussion in Section IIII Bl First of all, because of the 
(small) orbital eccentricity, the numerical flux oscillates around a "mean" value given by the PNQC estimate. The 
numerical flux starts deviating quite clearly from the 2PN and 3.5PN fluxes ~ 20 — 40M before the ISCO, and the 
agreement between numerical and analytical fluxes gets slightly better for larger mass ratio. A remarkable feature of 
the numerical flux is that it does not reduce to zero after the exponentially decaying ringdown phase. We believe this 
to be, at least in part 7 , an artifact of the memory effect discussed in Section Hi Al 



The 2.5PN expansion of the flux (not shown in the plots) has a physically unreasonable zero crossing ~ 20 — 40M before the radiation 
peak. The poor quality of the 2.5PN flux is in agreement with well-known results in the extreme mass ratio limit (see eg. Fig. 1 in lid ). 
Some contribution to the non-zero flux at late times may come from numerical noise. 
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FIG. 14: Total energy flux computed numerically (thick solid line), and substituting two different estimates of the orbital 
frequency (u>Dm and ui c ) into the 2PN and 3.5PN energy fluxes, i.e., keeping different numbers of terms in Eq. (|3.21|l . The 2PN 
and 3.5PN fluxes are so close they cannot be resolved "by eye" on the scale of this plot. The plots refer to runs D8. On the 
left we show results for q = 2.0, on the right for q — 3.0. 




FIG. 15: Dominant multipolar components of the energy flux computed numerically (thick solid lines), and substituting two 
different estimates of the orbital frequency (uiDm and ui c ) into the 2PN multipolar contribution to the energy flux, Eq. (|3.24[) . 
The plots refer to runs D8. On the left we show results for q = 2.0, on the right for q = 3.0. 



In Fig. [15] we look at the dominant multipoles contributing to the total energy flux: I = m = 2 and I = m = 3. 
The available analytical expansions of the mode amplitudes Mripi tTn are only 2.5PN accurate. Since the 2.5PN term 
in the I = m — 2 waveform depends on our particular choice of w, we dropped all terms of order higher than 2PN in 
Eq. (1X241) . 

By comparing Fig. [14] and Fig. [15] we can easily check that the I = m = 2 component carries almost half of the 
total energy flux (as long as we compare numerical results with numerical results, or consider the same truncation 
order in the PN approximation). The agreement between the 2PN estimate and the numerical flux is quite good for 
the dominant (I — m — 2) component. However, the 2PN approximation seems to systematically overerestimate the 
flux in I — m — 3. Given the present accuracy of our numerical code (and the poor convergence properties of the 
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total energy flux at 2.5PN) it is hard to tell whether this is an artifact of eccentricity in the simulations, or a genuine 
indicator of the convergence properties of the PN approximation for higher multipolar components. 

Whether we consider numerical results or the PNQC approximation, the relative contribution of higher-Z multipoles 
to the flux can be seen to increase with mass ratio. It also increases (for fixed mass ratio) as we get close to merger. 
To leading order in a PNQC expansion we can show that the ratio of the dominant multipolar components of the flux 
as a function of frequency (or alternatively, as a function of time to coalescence) is given by 
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FIG. 16: Relative contribution of different multipolar components to the integrated energy flux as a function of time. The plots 
refer to runs D8. On the left we show results for q — 2.0, on the right for q = 3.0. 



In Fig. [TB] we show the ratio of the integrated (numerical) energy flux in different multipolar components as a 
function of time. This plot confirms that the relative contribution of higher multipoles increases for large mass ratio, 
and (for given mass ratio) it increases as we get close to merger. 

Finally, in Fig. [17] we compare the numerical angular momentum flux with the PNQC prediction. The oscillations 
in the numerical flux seem to be a general feature: compare eg. Fig. 27 of BCP, where they are attributed in part 
to improper initial conditions in the time integrals required to obtain the flux from ^4. Our results confirm that, as 
remarked by BCP, Eq. (|3.25p seems to hold on average throughout the whole inspiral (possibly with larger deviations 
close to merger). In addition we point out an interesting correlation between the energy and angular momentum 
fluxes. In Fig. 1171 besides the angular momentum flux, we also plot the energy flux (multiplied by 50 for scale). 
The plots clearly show that oscillations in Fj have (roughly) the same period as oscillations in i*E • Perhaps this could 
be evidence that the observed oscillations are somehow related with the orbital eccentricity, minima and maxima 
corresponding to periastron and apastron. A detailed study of this correlation is beyond the scope of this paper. 
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FIG. 17: Total angular momentum flux computed numerically (solid line) and substituting two different estimates of the orbital 
frequency, Uom and u c , into Eq. (|3,25[l . We overplot the energy flux (multiplied by 50) for run D8 with q — 2.0, to show that 
zero-crossings of the angular momentum flux correspond (roughly) to local extrema of the energy flux. 
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IV. THE MERGER- RINGDOWN TRANSITION 



The goal of this Section is to study the ringdown phase, and to explore the properties of the final black hole formed 
after merger. We will compare different fitting methods to extract information from the ringdown waveforms. As 
discussed in [T3|, such a comparison can help us resolve real physical effects (such as, for example, time variations of 
the ringdown frequencies) from systematic parameter estimation errors due to the variance and bias of each particular 
fitting algorithm. In particular, here we consider two classes of fitting algorithms: the matrix pencil (MP) and 
Kumaresan- Tufts (KT) methods, which are modern variants of the so-called Prony linear-estimation algorithms for 
damped exponentials in noise [33l l34j ; and a standard non- linear least-squares technique, the Levenberg-Marquardt 
(LM) algorithm [Hj]. 

In [I7J we pointed out that Prony methods have a number of advantages with respect to standard non-linear least- 
squares techniques: (i) They do not require an initial guess of the fitting parameters; (ii) They provide us with a 
simple, efficient way to estimate QNM frequencies for the overtones, and even to estimate how many overtones are 
present in the signal; (iii) Statistical properties of Prony-based methods in the presence of noise (such as their variance 
and bias) are well studied and under control. When compared with the LM algorithm, Prony methods seems to have 
comparable variance but slightly smaller bias. 

In BCP, the real and imaginary parts of ipi iTn were fitted separately using standard non-linear least-squares methods. 
Prony-like methods allow us to fit the "full" , complex signal by a function of the form 



where <l>i' m ' n (j, M& n ) denotes a complex QNM frequency. In BCP the final black hole's mass and spin (j , Mfi n ) are 
taken as the independent fitting parameters, and the different QNM frequencies Coi mn are obtained, for given (j, Mg n ), 
either by using fitting relations or by interpolating numerical tables .36:] . 

BCP allow for general mode-mixing due to the expansion of spherical harmonics in terms of spheroidal harmonics. 
We will assume that each spherical (/ , m) mode is well described by a single (I , m) ringdown mode. Another difference 
is that BCP include overtones in the QNM expansion. Adding overtones provides a good fit of the strong-field phase 
by effectively increasing the number of fitting parameters (mode amplitudes Ai m n and phases 4>i mn of the overtones). 
This idea is perfectly consistent with QNM expansions in the context of linear black hole perturbation theory. An 
obvious drawback of the idea is that it assumes the validity of linear perturbation theory to extend the QNM fit 
before the peak of the radiation. Another potential problem is that, by using many fitting parameters, we can always 
get very good agreement with the numerical waveforms, but we do not necessarily get a better physical description of 
QNM excitation. For simplicity, in this paper we do not attempt to include overtones in the fit, but we only assess the 
accuracy of fits of the fundamental QNM. In [17] we have shown that the QNM frequency and damping time evolve 
quite rapidly right after merger. This evolution could be interpreted as a bias in the fitted frequencies induced by 
the omission of higher overtones; or, alternatively, it could mean that the mass and angular momentum of the newly 
formed, dynamical black hole spacetime really are evolving on timescales much smaller than the QNM timescales, 
producing an effective redshift in the QNM frequencies [33, [38[ . Issues such as the inclusion of overtones and the 
detailed study of nonlinearities will be addressed in the future. 



Independently of the chosen fitting method, there is some arbitrariness in choosing the time interval [to, tf] to 
perform the fit. A well-known problem with the merger-ringdown transition is that we do not know a priori when the 
ringdown starts P, [H, H^]- This problem is discussed at length in Section TlV Dl below. Ideally, the starting time for 
the fit to should be determined by a compromise between the following requirements: (i) to should be small enough 
to include the largest possible number of data points: in particular, we do not want to miss the large amplitude, 
strong-field part of the waveform after merger; (ii) to should be large enough that we do not include parts of the 
waveform which are not well described by a superposition of complex exponentials: the inclusion of inspiral and 
merger in the ringdown waveform would produce a bias in the QNM frequencies. 

A judicious choice of tf is also necessary. Usually we would like the time window to be as large as possible, but 
Fig.[T]and Fig. [5] clearly show that the low amplitude, late-time signal is usually dominated by numerical noise (mainly 
caused by reflection from the boundaries). This noise can reduce the quality of the fit, especially for the subdominant 
components with I > 2 and for large values of to- A practical criterion for the choice of tf is suggested by a look at 
Fig. [21 If the ringdown waveform were not affected by noise from boundary reflections, \Mr ipi <m \ should decay linearly 




(4.1) 



I'rn'n 



A. Choice of the fitting window 
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on the logarithmic scale of the plots 8 . At low signal amplitudes, we see boundary noise-induced wiggles superimposed 
to this linear decay. The first occurrence of these wiggles is a good indicator of the time t f at which numerical results 
cannot be trusted anymore. To test the robustness of fitting results to late-time numerical noise, while at the same 
time keeping the largest number of data points in the waveform, we decided to use two different "cutoff criteria" : 

1) "Relative" cutoff: remove from the waveforms all data for times t > tt — t le \, where t Ie \ is the time when 
the amplitude of each multipolar component \Mripi im \ becomes less than some factor ip cu toS times the peak 
amplitude \ipi ,m{t P eak)\ (values of i pea k for I = m = 2 are listed in Table [TTJ) : 



\Mr^, m (t iel )\ 

< ^cutoff- (4.2) 



I Mr V>;,m(ipeak)| 

2) "Absolute" cutoff: remove from the fit all data with t > tf = £ a bs, where £ a bs is the time at which the absolute 
value of the amplitude \Mrtpi tm \ < ^cutoff/10. 

The choice of the cutoff amplitude is somewhat arbitrary. We chose ^cutoff — 10~ 3 for low resolution, and V'cutoff = 
10 -4 for high resolution. 

For each chosen tf, we compare the different fitting routines as we let to vary in the range [t pe ak, tf]. By monitoring 
the convergence of the QNM frequencies to some "asymptotic" value as to — > oo, we can tell if the black hole settles 
down to a stationary Kerr state, or if, on the contrary, non-linearities and mode coupling are always present. Notice 
that as to grows the signal amplitude decreases exponentially, and we effectively reduce the signal-to-noise ratio (SNR) 
in our fitting window. Robust fitting methods should give reasonable results even for large values of t (that is, modest 
values of the SNR). 



B. From ringdown frequencies to black hole parameters 

In [ItJ we fitted the frequency ui and quality factor 9 Q of the I = m = 2 fundamental mode of the newly-formed 
black hole as a function of to- The results show that the QNM frequencies evolve quite rapidly in the first 10M — 20 M 
after merger: see in particular the bottom panels of Fig. 7 in [17] . where a rapid decrease of Q{to) is clearly visible 
for simulation times 240 < to/M < 260 M. Assuming linear perturbation theory to be valid, the real and imaginary 
parts of each QNM frequency are unique functions of the mass Mfi n and of the (dimensionless) angular momentum 
j = Jfl n /-M| n of the final black hole: say, M fila u>; TO „ = fi mn (j), M fin a; m „ = M &n /n mn = gi mn {j) HI]. The quality 
factor of the oscillations Qi mn , being dimensionless, must be a function of j only. A numerical calculation shows that 
for the dominant modes (I = m = 2, 3, 4) this function is monotonic and invertible (see eg. Fig. 5 in Ref. [36|). 
Therefore we can easily invert Qi m n{j) to compute j{to), either by using fitting relations or by interpolating QNM 
tables. 

The results of this inversion for the fundamental I = m = 2 mode of the black hole formed as a result of inspirals 
with q — 1.5 and q = 3.0 are shown in Fig. 1181 As the origin of the time axis we choose the time £ pca k at which 
the I = m = 2 amplitude has a maximum (see Table IP!)) . Solid lines refer to the "absolute" truncation criterion, 
and dashed lines to the "relative" truncation criterion (see Section TlV Ap . On the scale of these plots, different 
truncation criteria affect the estimated parameters only for low-resolution simulations and at relatively late starting 
times {to/M > 50), when the signal amplitude becomes comparable to numerical noise. Not surprisingly, there is 
remarkable agreement between KT and MP methods. The main difference when we use the non-linear least-squares 
LM method is a systematic time-shift in the angular momentum: the blue lines would be in excellent agreement with 
the prediction from Prony methods if shifted backwards in time by Ato ~ 2 — 3M. This time shift can easily be 
understood. In the non-linear least squares fit we are ignoring the imaginary part of the waveform. Since the real and 
imaginary parts of the waveforms are essentially time-shifted copies of each other, this produces a constant dephasing 
in the predicted physical parameters of the final black hole. 



With larger resolution and longer running times, eventually the exponential decay should turn into the well known power-law tail induced 
by backscattering of the radiation off the spacetime curvature |40j . In the simulations we consider, noise produced by boundary effects 
is large enough that this effect is not visible. 

We recall that the quality factor Q = \w/(2a)\, where a = X/t is the imaginary part of the QNM frequency, i.e. the inverse of the 
damping time. 
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FIG. 18: Estimate of the angular momentum from a fit of the I — m — 2 waveform using different methods. Top panels refer 
to a merger with q — 1.5, bottom panels to a merger with q = 3.0. Results on the left were obtained from low-resolution D7 
runs, and those on the right from high-resolution D7 runs. 



In the absence of numerical errors and mode coupling j(ta) should monotonically decrease, approaching a constant 
as to — > oo. Fig. 1181 clearly shows that this is not the case. All fitting routines consistently predict non-trivial time 
variations (roughly of order a percent) in j. Increasing the resolution reduces the amplitude of these variations, and 
produces a flattening of j(to) for 40 < to/M < 60. The angular momentum increase that can be seen for q = 1.5 and 
(to — tpcak) ^ 45M, and the oscillations in j for q = 3.0 in the same time range, are clearly artifacts of insufficient 
resolution. We tried to perform a Richardson extrapolation of the results assuming second-order and fourth-order 
convergence, to determine if angular momentum oscillations (which could be a sign of "new" physics) disappear in the 
limit of infinite resolution. Our results are shown in Fig.[TjJJ They are compatible with the possibility that oscillations 
disappear in the limit of infinite resolution, but more simulations and better control of the errors are required to reach 
a firm conclusion. 

Fig. [TS] clearly illustrates that resolution plays a role in the accuracy with which we can estimate black hole 
parameters from ringdown fits, especially at late times, when the signal is very weak and affected by numerical noise 
(likely caused by reflections off refinement and outer boundaries). Fortunately, changing the extraction radius does 
not affect the quality of the fits. We checked this by fitting the I = m = 2 and I — m = 4 modes for equal mass 
(q = 1), large separation (D10) binary mergers with different extraction radii r ex t = 30, 40 and 50. The functional 
form of j(to) is exactly the same at different extraction radii. Changing r cxt only produces a trivial shift of the time 
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FIG. 19: Richardson extrapolation ol the estimated angular momentum for I = m = 2 and q — 1.5, assuming second-order 
(thick line) and fourth-order (thin line) convergence. 



axis by Ato ~ Ar ext , due to the finite propagation speed of the waves. 

Estimates of the angular momentum as a function of to, obtained by fitting the dominant mode (I = m = 2) for 
different values of q, are shown in Fig. 1201 The angular momentum is constant within about ~ 1%, but the quality of 
the estimates rapidly degrades with mass ratio. Even with high-resolution runs, the estimated angular momenta have 
errors ~ 10% for 9 > 3. In the next Section we will show that improved estimates are possible if we cross- correlate 
information from different multipoles of the radiation, making use of the no-hair theorem of general relativity. 



o.(,o 



0.40, 



l=m=2 



10 



20 



30 



40 



50 



Ml 



0.70 



O.l.ll 



0.40, 



l=m=2, High resolution 



10 



20 



30 



40 



50 



(.0 



FIG. 20: Angular momentum estimated applying the MP method (thick lines) and the LM method (thin lines) to the I = m = 2 
waveforms. Lines from top to bottom refer to different mass ratios: q = 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0. 



C. Cross-correlating information from different multipoles to determine the black hole parameters 



We already pointed out that the quality factor of each QNM Qi mn , being dimensionless, must be a function of 
j only; and for the dominant modes (I = m = 2, 3, 4) this function is monotonically increasing, so we can easily 
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invert Qi m n{j) to compute j(io) by using fitting relations or by interpolating QNM tables 36]. If linear perturbation 
theory were an exact description of the final black hole's dynamics, the value of the angular momentum obtained from 
different QNMs - that is, from different values of (I, m, n) - should be the same for all modes and all values of to. In 
practice this is only approximately true. First of all, non-linear effects should be present close to merger, so that linear 
perturbation theory provides only an approximation to the "true" oscillation frequencies (if the definition of QNMs 
makes sense at all in the non-linear regime). Secondly, mode mixing induced by the use of spin- weighted spherical 
harmonics with some given given (J, , m) rather than spin- weighted spheroidal harmonics, will produce additional QNM 
frequencies 10 with different Z's and the same m. Finally, numerical error and the omission of overtones will inevitably 
produce some bias in the estimation of the frequencies, whatever fitting routine we use to extract them. 

All of these effects should be reasonably small, especially at late times, since linear perturbation theory can be 
expected to be a good approximation once the final black hole is "reasonably close" to a Kerr state (where "reasonably 
close" is here a loosely defined concept that can be made more precise, for example, through the use of quantities 
such as the "speciality index" S [Ilj]). Conversely, if we estimate angular momenta by fitting different multipolar 
components of the radiation, we can determine when perturbation theory is a good description of the system by 
looking for points (or intervals) in time when the angular momenta obtained from different fits agree with each other. 




FIG. 21: Consistency in the radiated angular momentum, as predicted by different multipoles, for a, q — 1.5, D7 run. Thick 
lines use the MP method, thin lines use the LM method. Solid lines are obtained by removing all points for which the absolute 
amplitude drops below 10 -4 , and dashed lines by removing all points for which the relative amplitude drops below 10 -3 the 
peak value. The choice of fitting method and truncation criterion has negligible influence on the results. Arrows mark the last 
point in time when the fit can still be considered reliable, and different multipoles agree with linear perturbation theory of Kerr 
black holes. 



In Fig. [21] we plot the angular momenta estimated by fitting the dominant multipolar components of the radiation 
emitted in a q = 1.5, D7 merger. Angular momenta from I = m = 2 and I = m = 3 are generally in good agreement, 
but they display oscillations around some mean value. The magnitude of the oscillations is larger for I — 3, and it also 
gets larger for coarse resolutions. However, there are discrete points in time when the angular momenta predicted by 
different multipolar components agree with each other. 

In Fig. [5Uwe use QNM fits of different multipoles to extract the final black hole mass Mg n . From Mfi n we can 
estimate the radiated energy as a function of to by computing (M — M^ n )/M. The plots provide a remarkable 
consistency check of the results in Fig. 1211 whenever results from numerical relativity are in agreement with linear 
black hole perturbation theory for the angular momentum, they are also in agreement for the radiated energy. In 
other words: when angular momenta from I = m = 2 and I = m = 3 agree, also the masses do. In our opinion this 
result is non-trivial, and it lends support to choosing this "perturbation theory time" (marked by arrows in the plots) 
as our best guess to estimate the final black hole's parameters. 



This mode mixing was actually observed by BCP: when fitting the I = 3, m = 2 waveform they also found the I = m = 2 QNM 
frequency. 
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FIG. 23: Angular momentum estimated applying the method to I = m = 2 waveforms (thick lines) and toi = m = 3(/ = ra = 4 
in the case q = 1.0) waveforms (thin lines). Hollow circles mark the "perturbation theory time" for each mass ratio (see text). 
From top to bottom, different linestyles refer to q = 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 and 4.0, respectively. 



In Fig. [221 we show the performance of the MP method in estimating angular momenta for different mass ratios: 
q = 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0. Increasing the resolution produces a flattening of all curves, the effect being 
more pronounced for large mass ratios. Remarkably, we find that the angular momenta and masses predicted from 
fitting different multipolar components agree at some well-specified time for all mass ratios. Now, in linearized 
theory different multipoles should be consistent with a single linearly perturbed Kerr black hole. This means that the 
predictions for M or j from one multipole should agree with the corresponding predictions from any other multipole: 
lines in Figs. [H] or [221 corresponding to different multipoles should lie on top of each other. Accordingly, we will take 
the latest (in time) of these points as our "best guess" to estimate the parameters of the final black hole, since by then 
the background dynamical spacetime is as close as possible to a stationary Kerr solution. In Table [V] we list the final 
angular momenta and radiated energies extracted from a QNM fit at this "optimal" time, comparing results against 
the corresponding estimates from wave extraction techniques. The values thus obtained from the two independent 
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TABLE V: Energy radiated -Eqnm/M and final dimensionless angular momentum j'qnm as computed by a QNM fit of high- 
resolution D7 runs. We use the "absolute" cutoff criterion for the fits. Results are presented using two fitting methods: MP 
and LM (the latter results are in parentheses). For ease of comparison, we also give the energy radiated as computed by wave 
extraction techniques. We denote by tpT the "perturbation theory time" when the I = m = 2 and I = m = 3 predictions for 
mass and angular momentum are in agreement. For comparison we also show estimates of the radiated energy Etot/M and of 
the final angular momentum j& n obtained by subtracting the radiated energy and angular momentum from the initial ADM 
mass and angular momentum. Uncertainties are estimated using Richardson extrapolation as described in Sec. [TTJ 
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methods agree within the error estimates, indicating that QNM fits facilitate estimates of the radiated energy and 
the final spin accurate to within a few percent or better. 

D. Criteria to determine the ringdown starting time 

There are important motivations to try and define the ringdown starting time and to isolate, in a non-ambiguous 
way, the energy radiated in the ringdown phase. For instance, from a detection-based point of view, the SNR of a 
ringdown signal scales with the square root of the energy in the signal [3ol l42l ] . To define the energy in ringdown 
waves we must somehow define the ringdown starting time. Being able to define the ringdown starting time is also 
important when comparing numerical simulations with PN estimates of the energy, angular and linear momentum. 
In fact, it has been suggested that the discrepancy between PN estimates and numerical results for black hole recoil 
is due to neglecting the ringdown in the former [3j. To check the validity of this statement we must, again, define the 
starting time of the ringdown phase. 

Unfortunately, early studies in quasinormal ringing have established that there is no such thing as "the" ringdown 
starting time (see eg. [39( and references therein). In fact, the waveform can never be exactly described as a pure 
superposition of damped sinusoids: it is always contaminated by noise or by other contributions (such as prompt 
response or tails). This is essentially a consequence of the incompleteness of QNMs. However, from a practical 
viewpoint the signal is indeed dominated by ringdown at some stage, and this is the reason why we can use ringdown 
waves to estimate black hole parameters [3a |. The time span of the ringdown phase can be defined in different ways, 
depending on context. In the followin g w e will discuss and implement three possible alternatives, two of which have 
already been proposed in the past (i~8l Il9||. 

1. A least-squares approach 

A natural way to determine the QNM content of a given signal would be to perform a non-linear fit of the data 
to an exponentially decaying sinusoid. Here the unknown parameters are usually found in a least-squares sense, by 
minimizing some functional of the form J2 t=u [ h (t) ~ h QNM (t, {A})] . In our specific case h would be the numerical 
data, sampled at instants t — ti, and /i^ NM (i, {A}) is the model waveform (an exponentially damped sinusoid). 

The model depends on a set of unknown parameters {A} over which the functional should be minimized. It is of 
course very tempting to treat the starting time as one of those parameters. This is a possible way to determine the 
ringdown starting time, and it served as the basis for the proposal in (l8| . There it was shown that the quality of a 



34 



QNM fit can be monitored by using some suitably defined norm. In particular, Ref. [l8| proposed to use 

f tf Hl m (t) - i/l^dt 
\\N\\(to) = Jt ° , * . ( 4 -3) 

where ^ fi l m has been defined in Eq. (|4.ip . Clearly ||iV|| — * when the fit is very close to the numerical waveform. 
The idea is that the norm should have a local minimum when the "trial" starting time tq tends to the "true" starting 
time, tq — * to- 




FIG. 24: Norm (|4.3[) as a function of the trial starting time for the dominant components (I, m) — (2, 2), (3, 3), (4, 4) in a 
merger with q = 2.0 (we consider a D7 run here). Thick lines are obtained by fitting the frequency each time we change To. Thin 
lines use the following fixed values for the QNM frequency: Mlor = 0.51677, Mm = 0.08586 for I = m = 2, Mu)r = 0.82210, 
Mm = 0.08571 for I = m = 3, Mljr = 1.12152, Mm = 0.08577 for I = m = 4. 



This idea works well for the classical perturbation theory problem of Gaussian pulses scattered off a Kerr background 
[HI], but unfortunately it does not provide a very clear answer when tested on binary black hole merger waveforms. 
The norm 1 1 JV| ] (t"o) f° r a binary with q = 2.0 is shown in Fig. [2U where it is computed in two slightly different ways. 
The simplest way treats the QNM frequencies as known: their values can be obtained once and for all by using Prony 
methods or non-linear fits [17] . and kept fixed as we change To- The second method achieves a marginal reduction of 
the norm by fitting for the QNM frequency at each starting time To- 

^From Fig. [M] we see that the norm has some of the desired properties. First of all, it grows as the quality of the 
QNM fit degrades: for example, it is larger for the subdominant (I, m) components. In addition the norm grows, as 
it should, when we try to extend the fit to encompass the merger region, i.e. when (tq — i pC ak) ^ 10. 

We find that the functional (|4.3[) has a minimum for most, but not all of the waveforms. Even when it does have 
a minimum (as in the case of Fig. I24[) this minimum is very broad. In addition the norm oscillates with a period 
which is basically the QNM period, and it has a series of local minima and maxima. The broad minimum and the 
oscillations make it very hard to locate the starting time. Of course, the functional f|4. 3|) is by no means the only 
possibility. We experimented with some alternative functional forms of the norm, but with no success. On a positive 
note, the situation seems to improve when the waveform is computed at large extraction radii [43[. 



2. Nollert's Energy Maximized Orthogonal Projection (EM OP) 



A physically motivated notion of ringdown starting time time was introduced by Nollert [l9j . He realized that the 
problems with defining the starting time arise immediately at the onset: QNMs are not complete and not orthogonal 
with respect to any inner product, so a quantification of the energy (and therefore of the starting time) going into 
each mode, using standard "basis expansion" methods, is difficult (if not impossible). The lack of orthogonality can 
be circumvented by formally defining an ortho gon al decomposition of the waveform into the contribution of one (or 
more) QNMs, and some orthogonal remainder [19] : 



h = h\\ + h± 



(4.4) 
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Here, h\\ and h± are the part of h parallel and perpendicular, respectively, to a given QNM or a finite number p of 
QNMs. We therefore write 



/in 



5> 



(i) h (i) 
"•QNM 



(4.5) 



where 



'QNM 



if t < t 

e~ W;t sin (ui r t + (j)) if t > to . 



(4.6) 



is the QNM, assumed to start at some time to- The decomposition is achieved using a standard orthogonal projection 



(h\\,h ± )=0, 

where the inner product, following arguments by Nollert [l9| . is defined in an energy-oriented way: 



One can show that the energy "parallel to the QNM component of the signal" is given by 
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(4.9) 



It is now meaningful to talk about (say) "the fraction of energy going into the first QNM" . This fraction obviously 
depends on the starting time to in Eq. (|4.6p . Nollert observes that the ratio of the energy "parallel to the QNM 
component" to the total energy in the signal, E\\/E to t, has a maximum as a function of to. We can define the 
ringdown starting time as the time to corresponding to this Energy Maximized Orthogonal Projection (EMOP). In 
other words, according to Nollert's criterion, the ringdown starting time to = £emop is chosen by looking for 11 
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h {l) * halt 



h {i) * h {i) dt 
' l QNM ,l QNM u ''' 
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(4.10) 



The previous integral is evaluated separately for each polarization component. To avoid memory effects, when we 
integrate iff 4 we fix the integration constant so that h — at the end of the simulation. We denote by -Eemop the 
maximized energy parallel to the QNM component of the signal: 

Eemop = E\\(t = tEMOp) ■ (4-11) 

Using Prony methods or non-linear fits [l7T | we first determine the QNM frequency and the damping time (for 
simplicity we consider a single QNM). Then we compute £emop and -Eemop by maximizing (|4.10[) over both to and 



TABLE VI: EMOP data for 1 — 2. Numbers separated by a comma correspond to the + and x polarizations, respectively. The 
fraction of the total energy in the I = 2 mode is about 42% for all mass ratios. We find that, independently of mass ratio, the 
value of £emop for a given polarization is generally at a fixed position relative to the maximum of the waveform's amplitude 
tpeak- We measure this relative difference by AtEMOP = £ P cak — £emop, which turns out to be roughly independent of q. 
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Another conceivable definition would not use the total energy in the waveform -Etot, but the energy in the waveform for t > to. It turns 
out that this quantity does not have a well-defined maximum. It is also possible to use a variable frequency in l|4.10ll . in which case one 
could possibly obtain a larger maximum. This method would be equivalent to matched filtering, which is discussed below. 
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FIG. 25: EMOP for I — 2 computed using run D7 (and high resolution). In the left panel we overplot e+ and the actual 
waveform, marking £emop by a vertical dashed line. In the right panel we show that results are quite insensitive to q: each line 
corresponds to a different mass ratio, and linestyles are the same as in Fig. [23] 



Our results for I — m = 2 and run D7 are presented in Table PVTl and Fig. [23 In the plots, e + x is the fraction of 
energy radiated at t > to in each of the two polarization components, normalized to the total energy radiated in the 
simulation, and computed for the value of the phase maximizing the EMOP. The first thing to notice is that there is 
a sharp maximum of the fractional energy going into ringdown. As seen from Fig. [23 ~ 42% of the total energy in 
the I = 2 merger waveform goes into ringdown. The results differ (very) slightly depending on the chosen polarization 
state. 

In Table PVTl we measure the ringdown starting time £emop relative to the peak in \Mrip22\, i.e., we compute 
AiEMOP = ipcak - ^emop- We see that AiEMOP is basically constant for all mass ratios and for all runs, corresponding 
to different initial separation of the binary. This is an important consistency test on the results. Notice also that 
^emop is located before the peak location. 

EMOP times for the two polarizations are displaced by about 3M for run D7. To define a unique ringdown starting 
time we take the average of both polarizations (in Table IVI1 an average over the two polarizations is denoted by 
angular brackets). Using this average starting time we can define an energy radiated in ringdown, also shown in Table 
IVII We find the following formula to be a good fit for the energy in the I = 2 mode: 



TABLE VII: EMOP data for / = 3. In this Table, by "peak" we mean the peak in the amplitude of the I = 3 mode. Numbers 
separated by a comma correspond to the + and x polarizations, respectively. The fraction of the total energy in the I = 3 
mode is about 44% for all mass ratios. 



q 


run b emop 

Etnt 


(tEMOp) 
M 


Af EMOP 
M 


(AfEMOp) 
M 


io 4 £: EMO p 

JVl 


(io 4 b E mop> 
M 


1.5 


D7 0.44,0.45 


233.0 


4.4, 6.4 


5.4 


3.5,3.9 


3.7 


2.0 


D7 0.45,0.45 


230.5 


7.4,5.4 


6.4 


7.7,6.9 


7.3 


2.5 


D7 0.44,0.45 


232.5 


8.2,6.2 


7.2 


9.6,8.6 


9.1 


3.0 


D7 0.44,0.45 


234.0 


8.4,6.4 


7.4 


10,9.1 


9.6 


3.5 


D7 0.43,0.45 


238.2 


4.8,7.3 


6.0 


7.8,9.2 


8.5 


4.0 


D7 0.43,0.45 


240.2 


4.8,7.3 


6.0 


7.5,8.8 


8.1 


2.0 


D8 0.44,0.45 


456.5 


7.1,5.1 


6.1 


7.6,6.8 


7.2 


3.0 


D8 0.45,0.44 


412.7 


4.8,7.3 


6.1 


8.8, 10 


9.4 



Results for I = 3 follow the same pattern (see Table \VH§ . The average (£emop) for I = 3 is located about 6M — 7M 
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after the average (^emop) f° r I — 2- The following formula provides a good fit for the energy in the I = 3 mode: 

ggM2£ = o,io4 , 1 = 3. (4.13) 

If we take ^emop for I = 2 as the fiducial ringdown starting time, we can compute the energy, angular and linear 
momentum radiated during the ringdown phase (as described by the EMOP). The results of this calculation are listed 
in Table QJ and fitting formulas are provided in Appendix [B] (see in particular Table IXIip . 



3. A detection-based approach: energy deposited in matched filters 

As we already stated QNMs do not form a complete set, so the signal will always comprise quasinormal ringing 
plus some other component (such as prompt response or tails). However, in most practical applications we are only 
interested in some "fairly good approximation" to the ringdown waveform. The notion of "fairly good" must be 
defined according to the specific context. 

A possible definition, based on theoretical considerations, was introduced in the previous Section. Here we propose 
an alternative, practical definition of the ringdown phase from a detection perspective. Detection of ringdown waves 
is likely to be achieved through matched filtering [361 142||. The technique works by cross-correlating the detector's 
output against a set of theoretical templates. It can be shown that the maximum SNR is achieved when the template 
is equal in form to the detector's output (hence the name matched filtering). Matched filtering is the method of choice 
to search for ringdown waves: it is quasi-optimal and inexpensive, in the sense that it achieves the maximum SNR 
with a relatively small number of templates or filters. 

Now, for the purpose of a matched filtering detection, the ringdown definition must be related to the use of ringdown 
templates. The relevant question is therefore: what is the maximum SNR attainable through the use of a filter which 
is a pure damped sinusoid? By definition, given the numerical waveform h(t), the SNR p is 

, — - m^m (Ws2 rM(/)w)+M/)fts(/) > (4 . 14) 

W,t v /(r({A},to)|T({A},t ) Jo S h (f) 

where the template T({\} , to) is 

T({\} ,t ) = l r nt ~ t0) + ^ ' * * - *° ' (4-15) 

^ if t < to • 

Sh(f) is the noise spectral density of the detector and {A} is a set of parameters characterizing the templates. The 
procedure is now simple: we "slide" this template backwards (starting at large to and decreasing it progressively) 
across the numerical waveforms, and determine the maximum of the convolution ^. 14[) . A good initial guess for the 
template parameters {A} = (tof ,lu^ ,4> t ) can be obtained with Prony methods [17| . 

As expected to will depend on the observer, i.e., on the detector being used, through the noise spectral density Sh(f)- 
In practice, however, the dependence on the detector is usually very weak, since in general the largest contribution to 
the convolution integral is near the resonant frequency u> r . Thus, for all practical purposes, the detectors behave as 
if the noise were white: the spectral density Sh{f) can be approximated as constant and moved out of the integral. 
This assumption also allows one to sidestep the computation of the Fourier transform of the waveforms: by ParsevaPs 
theorem, the frequency integral can be turned into a time integral. A more complete analysis, taking into account 
the full structure of the detector's noise, is in preparation. 

A possible notion of effective ringdown starting time tMF according to a matched filter, which is useful to make 
contact with previous SNR calculations [36| . can be given simply as follows 12 . Define the effective starting time tMF 
as the instant for which 



Nollert's "theoretical" definition, explained in the previous Section, is not too dissimilar from a "detection-oriented" definition. Indeed, 
expression (139) in |l9l can be interpreted as the fitting factor between actual waveforms and ringdown templates (for white noise). If 
we take the ringdown frequencies as unknown parameters and choose them to maximize the EMOP 114.1011 . the results we get are very 
close to the present matched-filtering criterion. 
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where p is computed from Eq. (|4.14[) . Notice that, in general, £mf does not coincide with the instant at which the 
convolution between the signal and the template has a maximum. By using Eq. (|4.16[) the SNR can be expressed 
in terms of energy in the actual signal. This is a common approach in engineering, introduced in the context of 
gravitational wave detection by Flanagan and Hughes [42j (see also [36j]). 



TABLE VIII: Estimated, polarization-averaged effective starting times tMF and energy radiated in ringdown from a matched- 
filter detection perspective, as functions of mass ratio. The listed energies should be taken as rough estimates, depending on 
the number of filters one is willing (and able) to use. We list also A£mf> which is the "effective" starting time as measure from 
the peak of the I = m = 2 waveform: A4mf = tpeak — tyiF- 

q (*mf/M) {At MF /M} {Emf/M~ 



1.0 


207 


27.0 


0.028 


1.5 


208 


28.4 


0.026 


2.0 


209 


26.4 


0.021 


2.5 


209 


29.1 


0.018 


3.0 


210 


29.8 


0.015 


3.5 


211 


31.5 


0.012 


4.0 


212 


32.8 


0.011 



The detection-based criterion, when applied to the merger waveforms considered in this paper, yields the results 
shown in Table I Villi From the above discussion, it is clear that the values we list for the energy radiated during 
ringdown are effective energies measured by the detector. These correspond to the values used in data analysis (see 
for instance [36|). From the Table we see that the effective energy radiated in ringdown for an equal mass merger 
is ~ 3%, in very good agreement with the "guesstimate" by Flanagan and Hughes [42j], which has often been used 
in the literature to compute SNRs and measurement errors. We also note that this value is much larger than the 
energy estimated by the EMOP, typically twice as large. This happens because the filter is looking for the maximum 
correlation, usually implying that the best-match parameters (uj r and u>i) will differ significantly from the true signal 
parameters. 

Also notice that different polarizations yield slightly different energies and starting times. For instance, for equal 
mass mergers, we get tMF ~ 205 and £mf ~ 208 for the plus and cross polarizations, respectively. If we average over 
polarization states, this yields an effective radiated energy of ~ 2.8%. 

We also point out that the amount of energy depends on the parameter space to be searched. In principle, the 
correlation (j4.14[) is to be maximized over all possible values of u> r , u>i. In practice this would lead to a very large 
number of filters, so we must choose reasonable cutoffs on the parameters. For instance, in black hole ringdown 
searches one looks for modes with a quality factor typically smaller than ~ 20. It may be possible to increase the 
SNR and the amount of effective energy in ringdown by enlarging the parameter search (this would also allow us to 
search for ringdown modes of other objects, such as neutron stars or boson stars). A discussion of these issues will 
be presented elsewhere. 

To conclude this Section, we point out that a fit of the total effective energy radiated in ringdown, according to a 
matched filtering criterion, is: 

^1 « 0.44^^- . (4.17) 
M (1 + 9) 4 



V. CONCLUSIONS AND OUTLOOK 



The present study of binary black hole waveforms is, in many ways, only preliminary. The following is a partial list 
of important open problems. 

Using "hybrid" waveforms in data analysis 

The present study explored the physical properties of numerical waveforms and their relation with analytical 
methods. Our focus has been on providing analytical insight into the structure of the waveforms. For this reason, we 
deliberately avoided problems at the interface between numerical relativity and data analysis (see Ref. [44[ for some 
steps in this direction). We strongly believe that an analytical understanding of the numerical simulations will be 
useful, or even necessary, to bridge the gap between the (daunting) numerical task of generating waveforms, and the 
injection of these waveforms into a data analysis pipeline. 

The PNQC approximation studied in this paper provides a concrete example. We showed that the physical content 
of any given simulation can be reproduced quite accurately by substituting the orbital frequency in the dominant 
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waveform amplitudes, Eqs. (|3 . T[) . These "hybrid" PNQC waveforms can be used to create simple but accurate 
templates, and to interpolate between numerical waveforms with different physical parameters. 

Despite the recent progress in numerical relativity simulations are still computationally expensive. Hybrid template 
families could be injected in LIGO, or used in connection with LISA simulators in future rounds of the Mock LISA Data 
Challenges [451 ]. Semianalytical waveforms may significantly reduce the number of simulations needed for detection 
and parameter estimation, and they should be particularly useful when spins are included in the model. 

Removing spurious eccentricity and including additional physical parameters 

Our study clearly shows that the simulations have some small, but non-negligible, eccentricity. The eccentricity 
shows up as a typical modulation of all physical quantities of interest: the punctures' orbital velocity (Fig. [6]), 
the binary's orbital frequency (Fig. [7]), the energy and angular momentum fluxes (Fig. [TH and Fig. [TT)) . and so on. 
Measur ing, this spurious eccentricity, and possibly removing it by fine-tuning initial data, is an important open problem 
[ToL I27I [23] . Incidentally, the study of truly eccentric binaries could be relevant for massive black hole binaries to be 
observed by LISA [2l|. 

In the present study we completely neglect spins. There is mounting evidence, based for example on recent studies 
of binary black hole recoil, that spins will have a dramatic effect on the inspiral-merger-ringdown transition. An 
extension of our study to spinning, precessing black hole binaries is urgently needed. 

Stitching numerical and analytical waveforms 

For reasons of space, we decided not to address the important problem of comparing the PN phase evolution with 
the numerical phase evolution. This problem is central to connect the early inspiral phase with the merger phase, and 
it is a topic of active investigation. Since numerical evolutions show signs of eccentricity, comparisons of the phase 
evolution may benefit from the inclusion of eccentricity in the PN models as well. 

Another active research field concerns the problem of "stitching" PN and numerical waveforms. For the purpose of 
this stitching, do we need the full PN waveforms, or does the restricted PN approximation (including PN corrections 
in the phase, but not in the amplitude) work well enough? Does the number of cycles to be simulated numerically 
depend on the mass ratio and other physical parameters (eg. the spins)? We plan to return to these problems in the 
future. 

Bridging the gap with black hole perturbation theory 

Computational resources and resolution limitations reduce the accuracy of numerical simulations for large mass 
ratios. Unfortunately, many astrophysical black hole binaries could have mass ratio q = 10 or larger (see eg. [2l[ 
and references therein) . It is important to determine the maximum value of q that should be simulated in numerical 
relativity, or equivalently, the smallest value of q for which black hole perturbation theory can be considered adequate 
for detection and/or parameter estimation. In Appendix [C] we collect some results that may be useful in this context. 

Astrophysics and gravitational wave detection 

The most interesting applications of our results should be in astrophysics and gravitational wave detection. For 
example, the multipolar analysis of the radiation performed in this paper can be used to determine the cosmological 
distance at which we can test the no-hair theorem with LISA, LIGO or Virgo. Future extensions of this analysis to 
spinning binaries could also predict the parameter range (mass ratio, spin magnitudes and directions) in which the 
recoil velocity is astrophysically relevant, and the probability for these regions of the parameter space to be populated 
in astrophysical scenarios. 
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APPENDIX A: MULTIPOLAR DECOMPOSITION OF THE POST-NEWTONIAN WAVEFORMS 

Here we list the spin- weighted spherical harmonic components of a PN expansion of the Weyl scalar For 1 = 2: 
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55j7-107,„,^ n o/o „ ,„,^ 2173 + 7483*7 - 2047r? 2 ,„,^ n4/ o 
1 H ^ (Mfi) 2/3 + 2?r(MO) - -(Mil) 4 ' 3 



42 



1512 



-107 + 347? / 112\ , _ , 
2l 71 + [ w + T" J 111 ' ( M °) 



Mr ^2 ,ie 



8 /tt £M . 
3V5^ (MO) 



43 509 cy 9 \ r ^ x4/ c, 

77 H 7i 2 (Mil) 4 ' 3 

126 126 ' 168 ' J ' 



1 + 



79 



2077 - 17 
28 



{Mil) 2 ' 3 + 



27r-i(l + lnl6) 



Mil 



(Ala) 



(Alb) 



For I = 3: 



27 Vr^ (MO)3 



l-(4 - 277)(Mfi) 2/3 + 



i, My - 6 In (3/2) 



Mil 



,'123 1838 887 ,\ ,,^. 4/ o 



32 /tt 



Mrikflj* = -J-r,(Mil) 



,10/3 



(1 _ 3??) _ 193 -725, + 365,3 (Mn)2/3 



- ( -32. + 27T- y(57r- Hi) ) Wl> 



1 - 



2(7/ + 4) 



(MO)2/3 + 5,-7(7+10^2) 



607 _ 136 247 2 , f ; ) 
198 _ ~99 V ~ 198 77 ' 



4/3 



(A2a) 



(A2b) 



(A2c) 



For I = 4: 



Mr*.** - g£(M»f/- {(1 - 3,) - 1779 - 63 ^ + 



+ 



42 

4tt I — - 81n(2) 



162 FT SM 
-Vl4^^ ' 



/ 1 193 \ 
//( 127r-i^— -24In(2)J 



(1 - 277) 



39 1267 131 2 , , ,, 

TT + T32- ?? -^ H M °) 



2/3 



00 

Mr^ 4 , 2 e^ = ^-V^v(Mil) 10/3 
Do 



21 



(1 - 377) - 



1311 - 402577 + 285^ 
330 



(Mil) 2 / 3 



+ ( — (l- 477)i + 27r(l-3r?) ] Mil 



(1 - 277) + 



Mr^e* = ^V^iMiir/ 3 



101 337 88 . , , 
^3- + 14"^ -33^ ) (M0) 



2/3 



(A3a) 
(A3b) 

(A3c) 
(A3d) 



41 



For I = 5: 
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Our is the same as the inclination angle 1 in Blanchet et al. [l3| . We recall that 
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and therefore the phase <p is defined up to an additive term mc, with c a constant factor. By fixing the constant to be 
c = 7r/2 we recover, in the limit 77 — > 0, Poisson's results from perturbation theory. If we include radiation reaction 
terms in the waveform by using Eq. (27) of we find vj — —24. If we neglect those terms by redefining the phase 
and using their Eq. (32), then w = —8/7. 

The spin- weighted spherical harmonic components of h + ,h x can be obtained from the corresponding components 
of * 4 , Eqs. l(Ala| - ([A"6d"]) . as 



( h +- ih x)l,m = -^2Q2^l,m- 



(A8) 



The resulting expressions do include the logarithmic corrections to the phase. They are valid up to 2.5PN, with the 
only exception of the I — m = 2 component, which is given in Eq. 
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APPENDIX B: ESTIMATES OF THE POST-PLUNGE ENERGY, ANGULAR MOMENTUM AND 

LINEAR MOMENTUM 

In this Section we give estimates of the energy, angular momentum and linear momentum radiated in the last phases 
of a binary black hole inspiral. 

TABLE IX: Energy, angular momentum and linear momentum emitted after the estimated time of CAH formation, as listed in 
Table HT1 The CAH formation time is measured relative to the peak of the I = m = 2 waveform: AtcAH = (ipcak — £cah) /M. 



q 


run 


Mcah/M Ecah/M Jcah/M" 


EcakM/ Jcah 


lO^P* ,CAH 


/M 10 4 P a ,CAH 


/M 10 4 P C ah/M 


1.0 


D7 


19.0 


0.0259 


0.1261 


4.869 











1.5 


D7 


18.2 


0.0232 


0.1193 


5.142 


-2.56 


-0.28 


2.58 


2.0 


D7 


17.9 


0.0193 


0.1004 


5.202 


-3.98 


-0.15 


3.98 


2.5 


D7 


17.1 


0.0156 


0.0852 


5.462 


-5.01 


-0.20 


5.01 


3.0 


D7 


16.6 


0.0128 


0.0716 


5.594 


-5.38 


-0.53 


5.41 


3.5 


D7 


15.7 


0.0105 


0.0580 


5.524 


-5.62 


-1.34 


5.78 


4.0 


D7 


14.3 


0.0086 


0.0442 


5.139 


-5.78 


-2.30 


6.22 


2.0 


D8 


17.7 


0.0192 


0.1055 


5.495 


4.15 


-0.23 


4.16 


3.0 


D8 


15.2 


0.0125 


0.0493 


3.944 


0.28 


-5.90 


5.91 



Table IIXI lists the energy, angular momentum and linear momentum radiated after the time of CAH formation, 
estimated as the point where the ratio of the radial and tangential puncture velocities v r /vt > 0.3 (see Section Hill) . 

Given the uncertainties in our estimate of the CAH formation we also consider another useful (if somewhat conven- 
tional) indicator of the regime of validity of PN expansions: the Innermost Stable Circular Orbit (ISCO). The ISCO 
is defined by the condition that the energy of the two-body system £ , which is a function of the orbital frequency f2, 
has a minimum: d£/dil = 0. Since £ is only known as a PN series in f2, the location of the ISCO depends on the PN 
order [321 ]. 

TABLE X: ISCO data using the 2PN Taylor expansion of the energy. The ISCO time is measured relative to the peak of the 
I = m = 2 waveform: Atisco = (i pca k — iisco) /M. 
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TABLE XI: ISCO data using the 3PN Taylor expansion of the energy. The ISCO time is measured relative to the peak of the 
I = m = 2 waveform: Atisco = (i pca k - iisco) /M. 
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0.0645 


0.191 


-2.70 


-2.78 


3.88 


4.0 


D7 


0.107 


28.1 


0.0106 


0.0568 


0.187 


-2.28 


-3.07 


3.82 


2.0 


D8 


0.120 


23.6 


0.0210 


0.0940 


0.223 


+2.75 


+1.26 


3.02 


3.0 


D8 


0.112 


25.8 


0.0146 


0.0872 


0.167 


-3.10 


-1.86 


3.62 



In Tables |X] and IXII we list the orbital frequency at the ISCO Mflisco computed by including terms in the energy 
function up to 2PN and 3PN, respectively. Notice that 3PN corrections lower the ISCO frequency for all mass ratios, 
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the reduction being larger for larger mass ratios 13 . We also list the time location of the ISCO (relative to the peak 
in the amplitude of the I = m = 2 mode). We identified this time location as the instant when the ISCO frequency 
equals the orbital frequency from our simulations, as estimated from the gravitational wave emission of the dominant 
multipolar component I = m = 2 (see Section UlI Bp : MOisco = Mujd2- The 3PN ISCO "absolute" location in terms 
of the total simulation time is also shown in Table [IIJ and it should be compared with the C AH formation estimates 
in the same Table. As q — ► 1 the CAH formation time is very close to the 3PN ISCO time. For large mass ratios, 
when one of the holes is very small, the difference is larger, as expected on physical grounds. 

Tables IXl and IXj also list the energy, angular momentum and linear momentum emitted after the ISCO, with the 
ISCO location estimated by PN methods. While the energy emitted is a robust quantity, with very weak dependence 
on the PN order, angular and linear momenta are very sensitive to a variation of the PN order from 2PN to 3PN 
(i.e., they are very sensitive to small variations in the starting time of the integration). The reason for this behavior 
is apparent from an inspection of Figs. [T4l and 1 171 While the energy flux is a smooth function, even in the strong-field 
region, the angular momentum flux is a strongly oscillating function of time. 

The functional dependence of energy and angular momentum on mass ratio q can be inferred by combining the 
multipolar decomposition of the PN expansion (Appendix [A| with Eq. (I3.24|) . We find that good fits to the total 
angular momentum, energy, and multipolar energy distribution in the dominant modes for times t > to — i 3 |co are: 

Jisco/M 2 | 4>to =j tot -^L, (Bl) 



q 2 q 2 q 2 {q ~ I) 2 

E lSG o/M\ t>to = etot ^ ~ q y , £isco,2/M| t>to = e 2 —-— , E ISCO ,s/M\ t>to = £3 ^ - ^ 6 

with the fitting coefficients listed in Table IXlTl 



TABLE XII: Fitting coefficients for the energy and angular momentum emitted after the ISCO and EMOP times. 



to 


jtot 


etot 


£2 £3 


3PN ISCO 


2.029 


0.421 


0.397 0.168 


EMOP 


1.173 


0.295 


0.271 0.104 



In the same Table, for comparison, we also list the corresponding coefficients for t > to = £emop • The ringdown 
and plunge phase are strongly related with each other, and the numbers are roughly proportional. According to the 
EMOP criterion ringdown always starts after the ISCO. Therefore the post-EMOP radiation of energy, angular and 
linear momentum is always smaller than the corresponding radiation after the ISCO. 

The PN expansion breaks down after the ISCO. An estimate of the linear momentum emitted after the ISCO, 
within PN theory, was obtained in [2O] by integrating the PN linear momentum flux along a plunge geodesic of the 
Schwarzschild metric. In [2fJ, the integration is performed all the way from the ISCO (r ~ 6M) to the Schwarzschild 
horizon (r ~ 2M). The energy and angular momentum radiated after the ISCO can be computed using the same 
method, and they were kindly provided to us by Clifford Will [47| . 

Results of a 2PN estimate of the energy and angular momentum radiated after the ISCO are shown in Fig. [21] for 
different mass ratios, along with different estimates of the corresponding quantities from our numerical simulations. In 
particular, we show numerical estimates of the energy and angular momentum radiated after the 2PN and 3PN ISCO, 
and after the CAH formation, as functions of the mass ratio. To check the robustness of our results against initial 
conditions we also considered two runs starting at larger initial separations (namely, D8 simulations with q = 2.0 and 
9 = 3.0). 

Some comments are in order. The radiated energy from the simple PN estimate is in surprisingly good agreement 
with numerical results, the agreement getting better as we increase the PN order used to estimate the ISCO location. 
The agreement is particularly good when we consider radiation emitted after the 3PN estimate of the ISCO and 
relatively large mass ratios. The agreement in the radiated angular momenta is much worse. This seems to be a 
general feature when comparing PN estimates against numerical simulations. For example, Fig. 2 and 3 in [27| show 
that the eccentricity required to match PN predictions for a binary's angular momentum against numerical calculations 



For comparison, the ISCO for point particles is at ro = 6M, or equivalently at an orbital frequency of Mfl = 6 3 / 2 ~ 0.068. Corrections 
to this point particle limit were worked out by Clark and Eardley l46ll . yielding a simple analytical estimate for large but finite mass 

ratios: r /M = 6q/(l + q) , Mfi = r~ 3/2 . For q = 4 this yields MS! ~ 0.095, not too far from the 3PN estimate of 0.107. 
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FIG. 26: Energy and angular momentum radiated in the plunge using the 2PN and 3PN definitions of the ISCO are compared 
against the simple estimate by Blanchet et al. 20] (BQW in the legend). All estimates were computed using the D7 runs 
(except for the inverse triangles, which refer to D8 runs). 



in quasiequilibrium is significantly larger than the eccentricity required to match the corresponding energies. In the 
present case, the disagreement is partially affected by the strongly oscillating functional dependence of the angular 
momentum flux (see eg. Fig. ITT)). This is confirmed by the fact that a relatively small change in the initial separation 
(using D8 runs instead of D7 runs) produces a significant change in the numerical estimate of the angular momentum 
radiated after plunge. Given the large uncertainties associated with both numerical and analytical estimates, we 
cannot draw reliable conclusions from the observed disagreement. 

It may be tempting to attribute the observed differences in the angular momentum to the fact that the PN estimates 
of [20| neglect the ringdown phase. We can naively try to correct for this effect by integrating the energy and angular 
momentum fluxes from the ISCO up to the CAR only. However, this will result in serious disagreement for both the 
radiated energy and angular momentum: they turn out to be extremely small, especially for mass ratio q —> 1 (in this 
limit the CAH formation time and the ISCO are very close, see Table ITT)) . Thus, ringdown alone cannot explain the 
disagreement. A more detailed analysis, possibly combining the PN approach and the close-limit approximation, is 
necessary. 

An interesting possibility is that the agreement between PN estimates and numerical results could improve if the 
PN integration is truncated at the light ring, instead of integrating all the way to the horizon (as originally done in 
[20|). The physical argument for truncating at the light ring is that most of the radiation emitted after r ~ 3M would 
be filtered by the potential barrier surrounding the black hole, and that this potential barrier (for Schwarzschild black 
holes) has a peak at the light ring (23[. Fig. [27J shows that truncating at the light ring sensibly improves the estimate of 
the linear momentum radiated after plunge, correspondingly improving the estimate of the total kick velocity. Given 
the uncertainties involved in the extrapolation, this may be little more than a coincidence. In any case this problem 
is worth investigation, given the potential astrophysical relevance of recoil velocities. 



APPENDIX C: MULTIPOLAR DISTRIBUTION OF RADIATION FOR EXTREME MASS RATIOS 

In this Appendix we collect the main results for the energy, angular momentum and linear momentum radiated 
by particles falling into (rotating or non-rotating) black holes. Our purpose is to provide a quick reference for 
the extreme-mass ratio limit of numerical relativity simulations, to be compared with present and future numerical 
relativity calculations of binaries with large mass ratio (and possibly spin). Research in this direction is already under 
way: for example, evolutions of large mass ratio binaries using finite element methods can be found in 48]. 
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FIG. 27: The kick velocity accumulated after plunge using the 2PN and 3PN definitions of the ISCO is compared against the 
corresponding estimates by Blanchet et al. [2fj] (BQW in the legend). In the two BQW estimates the integration is truncated 
at the horizon or at the light ring, respectively. 



1. The energy radiated by plunging particles: non-rotating black holes 

The first investigations of particles plunging into black holes began with Zerilli [49( , who laid down the perturbation 
formalism to analyze gravitational radiation from a point-like particle with mass m p around a Schwarzschild black 
hole with mass M ^> m p . His analysis was completed by Davis, Ruffini, Press and Price 50] (hereafter referred to 
as DRPP), who numerically computed the gravitational radiation generated when a small particle at rest falls from 
infinity into a Schwarzschild black hole. DRPP found that the total energy emitted in the process (in geometrized 
units) is given by 

E tot = 0.0104(m*/M) . (CI) 

Detweiler and Szedenits [5l[ and Oohara and Nakamura [52} generalized DRPP's results to particles plunging 
into a Schwarzschild black hole with non-zero orbital angular momentum. In the perturbation framework under 
consideration, the particle's trajectory as it plunges down the hole, with zero velocity at infinity, is described by 

d or /I dt 1 d( t> L Z (n ~, 

9 = " /2 ' d7=l-2MA' = (C2) 



dr 

Tr = ± 



L 



r 



1/2 



1 - (1 - 2M/r) 1 + -f . (C3) 



The particle has an orbital angular momentum J p — m p L z . When L z = the particle falls straight into the black 
hole. For L z between zero and AM, the particle spirals a finite number of times around the hole before crossing the 
event horizon. For L z — AM, the particle spirals an infinite number of times around the marginally bound circular 
orbit at r = AM. For L z > AM the particle never falls into the black hole, so we discard this case. 

In Table [XTTTI we show results from 52] for the total energy radiated in the first three multipoles (I — 2, 3, 4) as a 
function of L z . From the Table it is apparent that the total energy output grows with L z , and so does the energy in 
each multipole I. The relative energy output in each mode behaves somewhat differently: as L z grows, the percentile 
energy going into I = 2 decreases. The opposite happens for all other modes. 

For L z = AM the total energy is obviously infinite: the particle spirals an infinite number of times around r = AM 
and therefore radiates for an infinite time. The energy radiated as a function of I is well approximated by a simple 
function: Ei ~ ae~ bl m p /M, where the coefficients a, b (which are functions of L z ) are listed in Table [XlTTl The total 
energy radiated is well approximated by E to t ~ ae~ 2b (l — e~ b )~ 1 rrip/M (52|. Remarkably, the relative contribution 
of the I = 3 mode is always larger than 10%, and that of the I = A mode is always larger than 1%. As usual the I = 2 
mode dominates, with a relative contribution always larger than ~ 50%. 

Also shown in Table [XlTTl is the number of spirals the particle completes before entering the horizon. This number is 
useful for two reasons. The first reason is that, if the particle falls with angular momentum very close to the marginal 
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TABLE XIII: Energy radiated in each of the three lowest multipoles for a particle with mass m p and angular momentum m p L z 
falling from infinity into a Schwarzschild black hole (from [52]). We show the percentage radiated in each mode relative to 
the total energy radiated (as extrapolated from the data, which typically yields an error of less than 5%). We also show the 
number N = [A0|/(2-7r) of "laps" the particle performs before plunging. The coefficients a and b are defined in the text. 



L z /M 


a 


b 


N 


-Etot 


E 2 




% 


E 3 




% 


E 4 




% 





0.44 


2 





0.010 


9.1 x 10 


-3 


88 


1.1 x 10" 


-s 


10 


1.5 x 10 


-4 


1.4 


1 


0.21 


1.5 


0.15 


0.013 


1 x 10" 


2 


78 


2.3 x 10" 


■3 


18 


5 x 10" 


4 


4 


2 


0.22 


1.1 


0.32 


0.036 


2.4 x 10 


-2 


67 


8.1 x 10" 


■3 


22 


2.7 x 10 


-3 


7.5 


3 


0.36 


0.86 


0.55 


0.112 


6.4 x 10 


-2 


57 2.7 x 10" 


-2 


24 


1.2 x 10 


-2 


10.7 


3.5 


0.53 


0.76 


0.75 


0.218 


1.2 x 10 


-1 


55 


5.4 x 10" 


-2 


25 


2.5 x 10 


-2 


11.5 


3.9 


1.0 


0.71 


1.15 


0.485 


2.4 x 10 


-1 


49 


1.2 x 10" 


-1 


25 


5.8 x 10 


-2 


11.9 



value AM, it will complete many revolutions and radiate a huge amount of radiation. In this case the perturbation 
expansion would no longer be valid, and therefore we must make sure that N is not much larger than one. The second 
reason is that N gives us an estimate of how much of the output energy is due to the actual, almost radial plunge 
motion (see eg. Fig. [5]), and how much of it comes from the particle circling around the black hole. 



2. The energy radiated by plunging particles: rotating black holes 

The standard formalism for small perturbations of Kerr black holes was formulated by Teukolsky [53j . The equations 
decouple and separate, reducing to two coupled ordinary differential equations with a source term. In the case of 
gravitational waves emitted by particles plunging into the hole the source term diverges at the boundaries, so this is not 
the most convenient formalism (but see [51[ for a way to get around these difficulties). Using the alternative formalism 
developed by Sasaki and Nakamura 54} , a series of papers by Nakamura and co-workers (see eg. (EH, [H, [5(1 [57[ and 
references therein) examined the gravitational radiation emitted by point particles moving in the vicinities of a Kerr 
black hole. 



TABLE XIV: Energy radiated in each of the three lowest multipoles for a particle with zero angular momentum falling from 
infinity into a Kerr black hole along the equator, as a function of j = J/M 2 . Taken from Fig. 3 in [Hrjl ]. 
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% 
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9.1 x 10" 
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-3 
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-4 


1.4 
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■2 


1.5 x 10" 


-2 


83 


2.2 x 10" 


■3 


12 


3.9 x 10" 


-4 


2.2 


0.85 


2.3 x 10" 


■2 


1.9 x 10" 


-2 


83 


3.4 x 10" 


■3 


15 


7.3 x 10 


-4 


3.2 


0.99 


4.7 x 10" 


■2 


3.3 x 10" 


-2 


70 


9.6 x 10" 


■3 


20 


2.7 x 10" 


-3 


5.7 



TABLE XV: Energy radiated in each of the three lowest multipoles for a particle with angular momentum m p L z falling from 
infinity into a Kerr black hole with j — 0.85, along the equator. Taken from Fig. 5 in [55| . 



L z /M 
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-l 
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-2 
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-2 


23 
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3 
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0.65 
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-2 


0.06 
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-2 


73 
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3 


18 
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3 
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0.0 
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-2 
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-2 
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■3 
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-4 


3.2 


-0.8 
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-3 
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10 
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-4 


1.4 
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-2 


0.38 
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-2 


81 
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■3 


12 
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-4 


2.3 


-4.5 
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-12 
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-2 


71 
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3 
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The results for the total energy, as well as the energy radiated in each mode, are summarized in Tables IXI VI and IXVI 
These Tables refer to particles falling along the equator of the black hole. Starting with particles falling from infinity 
with zero orbital angular momentum, Table IXTVl shows a familiar pattern. The total energy increases with increasing 
j. For near-extremal black holes (j = 0.99) the total energy radiated is almost five times the Schwarzschild value. 
Again, the total energy going into each multipole I increases. As previously pointed out, the relative contribution of 
each mode increases with j for I > 2, but it decreases for I = 2. The I = 3 mode always contributes more than 10% 
of the total energy, and the I = 4 mode always contributes more than 1%. 



47 



Table IXVl shows results for a j = 0.85 black hole, for several values of orbital angular momentum L z . The total 
energy varies by two orders of magnitude between L z — 2.6M and L z = — 0.8M . The energy emitted is larger than 
10% for I = 3, and larger than 1% for I — 4. We also list the number of spirals before plunge. This number is always 
smaller than unity, so the results can in principle be interpreted as a plunging motion and applied to the merger 
phase. 

Not shown here is the contribution of different m's to the total energy. For large black hole rotation and large 
positive values of L z most of the energy goes into I — m modes. Negative-m modes emit a negligible amount of 
radiation in this regime. For negative L z the situation is different: all modes seem to be excited to comparable 
amplitudes. See [H[ for more details. 



3. Linear and angular momentum radiated by plunging particles 



Computing the linear momentum carried by gravitational waves is of great astrophysical importance. Coalescing 
binary black hole systems may abound in galactic disks and in the centers of galactic nuclei. Due to the emission 
of gravitational radiation the final black hole receives a "kick," i.e., it acquires a non-zero recoil velocity because of 
momentum conservation. Depending on the momentum emitted, the recoil velocity may be large enough to release 
black holes from the host galaxy. If so, gravitational radiation effects will have considerable observable consequences 
for astrophysics and cosmology, such as the depletion of black holes from host galaxies, the disruption of active galactic 
core energetics, and the ejection of black holes and stellar material into the intergalactic medium. In the following we 
briefly review some perturbative calculations of the recoil velocity. 



a. Linear momentum 



In a Schwarzschild background, the linear momentum carried by gravitational waves when a zero angular mo- 
mentum point particle falls into the black hole is |AP| = 8.73 x 10 _4 m 2 /M [57]. This leads to a recoil velocity 
v ~ 2.63(10m p /M) 2 km/s. For the general case of a particle plunging with non-zero angular momentum, the linear 
momentum carried by gravitational waves is well approximated by 



|AP| = 9x 10~ 6 



.L z 
> — 
M 



10 



m 2 p /M, 0<^<3.4, 



L 2 



4.5 x lO- 2 m 2 JM , 3.4 < — < 4 



This leads to a recoil velocity of 
v ~ 2.7 



1£, 
2 M 



130 (lOwp/M) 2 km/s, 



(10m p /M) 2 km/s, 0< 



3.4 < —j- < 4 . 

M 



Lz 
M 



< 3.4. 



(C4) 
(C5) 

(C6) 
(C7) 



These results are not very sensitive to the rotation of the black hole (see eg. Fig. 6 in [55|]). A hint to extrapolate 
recoil velocities to mass ratios close to unity comes from the Newtonian result: replace /J./M by f(q) (58j . where 



/(?) = Q 



The function f(q) has two extrema at 

3- V5 

q = 



0.38, / 



2 g ~ 1 



25v / 5 



= /min(g) 



2.62, / 



V5 N 



25V5 



/max(O') 



(C8) 

(C9) 
(C10) 
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We then get 



2 (L z \ 1 L z 
5 V M ) + 2M + 



Akm/s. 0<^<3.4, 



232 ^ km/s. 



3-4 < < 4. 



(Cll) 
(C12) 



6. Angular momentum 



Radiation of angular momentum demands either a non-zero orbital angular momentum of the particle, or a non-zero 
angular momentum of the black hole. For Schwarzschild black holes, Oohara and Nakamura found that, even 
though both the radiated energy .Etot and angular momentum AJ diverge in the limit L z /M — > 4, their ratio is to a 
good approximation independent of L z : 



M E t , 



AJ 



0.15, for \L Z /M\>1 



(C13) 



For rotating black holes, the angular momentum radiated depends on the relative sign between the rotation of the 
black hole and L z . We denote by x the ratio of radiated energy to radiated angular momentum: 



L z 



ME tc 
AJ 



(C14) 



This quantity is listed in Table IXVll (after Table I in [55j) for some values of L z and j. 



TABLE XVI: The factor x(j, L z ) = M E tot /AJ for some values of L z and j (from Table I in |Sa|). 
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1 


0.15 


1.5 


0.19 


1.3 


0.21 


1 


0.30 


2 


0.15 


2.25 


0.20 


1.95 


0.22 


1.5 


0.30 


3 


0.15 


3.0 


0.22 


2.6 


0.24 


2 


0.34 



4. Perturbation theory as a guide to numerical results for comparable mass ratios 

One of the most prominent features borne out of binary black hole simulations seems to be the absence of strong 
non-linearities: the potential barrier close to the black hole horizon acts as a very effective cloak, filtering out many 
non-linear features of the dynamics 59]. For this reason results from perturbation theory (i.e., ij > 1) can usually 
be extrapolated to the equal mass ratio case, yielding very good agreement with full blown numerical simulations. 
To quote Smarr, "the agreement is so remarkable that something deep must be at work" [6(j. This was also found 
to be the case for the head-on collision of two black holes. The perturbation theory result E tot = 0.0104(Mf /M2), 
with Mi < M 2 Hi] can be compared to the full numerical result for q = 1: E tot w 0.0013M = 0.001 x 16ry 2 M = 
0.0104 x (2M)ry 2 (see [Hj]). Simple scaling arguments applied to perturbative results work surprisingly well. 

On this basis, a natural conjecture is that particles plunging with large but sub-critical orbital angular momentum 
should describe reasonably well the final stages of a binary black hole inspiral [5l|. Indeed, some of the results 
discussed in the main text suggest that the merger of two equal mass black holes can be described by extrapolating 
results from perturbation theory. For instance, the final spin of the black hole can be predicted by using the small 
mass ratio approximation: see Eq. (|3.20[) and the related discussion. 

For inspiralling binaries evolving through quasi-circular orbits, the plunge ("merger") happens when the smaller 
body crosses the ISCO, even though this notion is not well defined for comparable-mass bodies. We argue that 
the merger phase should be reasonably well described by a particle plunging with an orbital angular momentum L z 
only slightly smaller than the marginal value L z = 4M (corresponding to a radius equal to AM). In this case the 
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trajectory resembles that of a particle going through the merger phase: a quasi-circular orbit followed by a plunge. It 
is important to specify how close L z should be to the marginal value. We look for orbits that complete barely less than 
one lap before plunging. This guarantees that the energy output is due only to a plunge trajectory; it also guarantees 
that the orbit was quasi-circular before the plunge. We therefore argue that, as seen from the Table, L z ~ 3.9M is a 
near-optimal value. An obvious objection is that the ISCO for point particles is not at r = 4M (the location of the 
marginally bound orbit), but rather at r = 6M. Fortunately, for more massive bodies an approximate ISCO can be 
defined, and it is usually closer to r — AM. We can thus try to extend these perturbative results to understand the 
merger phase. 

Extrapolating the point particle results presented in the previous subsections we get, for non-rotating holes: 

Tf = °- 485 7t4^4 > ^rr = - 15 ' Recoil = 232-/^ km -s" 1 . (C15) 

M (1 + 9) AJ JWmax 

For instance, the value 0.485g 2 /(l + q) 4 for the energy is found extrapolating the last row in Table IXiTlI by the 
substitution m^/M 2 — > q 2 /(l + q) 4 . From fits of the corresponding numerical quantities at the 3PN ISCO (which 
should work as a good reference point, for lack of a better guess) we get: 

-ElSCO n Q 2 M Eisco n1fl ISCO ion /(#) l -1 frucl 

M (1 + g) 4 A Jisco /(<?)max 

The two sets of values are reasonably consistent. The largest disagreement refers to the linear momentum radiated, 
and therefore to the recoil velocity of the final hole. 

APPENDIX D: POLARIZATION OF THE WAVEFORMS 

Getting information about the polarization content of a waveform is simple in the presence of a monochromatic 
wave. In more general settings, making statements about the polarization state is easier in Fourier space. Methods 
to compute the polarization content of a waveform solely from a time-domain analysis were presented in [6l| . Here 
we shall adapt these techniques for the case at hand. 

From the (real) polarization components h+(t) , h x it) 14 we can define the so-called analytic signal H+ , H x in the 
following way: 

H+ = h+(t) + iH+{t) , (Dl) 
H x = h x (t)+iH x (t). (D2) 

The imaginary part of the analytic signal is the Hilbert transform Ti.it) of the signal h(t), defined as 

nit) = - f +Ca ^dr, (D3) 

where the integral is taken as the Cauchy principal value. For reference we note that the Hilbert transform of smt is 
— cost and the transform of cost is sini. 

From the analytic signal we define a covariance matrix C as 

( H + H* + H + H* \ , s 

{H X H* + H x H* x J ' {UV 

where an asterisk stands for complex conjugation. We note that the covariance matrix is Hermitian, and thus its 
eigenvalues Ao , Ai are real and positive. Without loss of generality, we assume Ao > Ai. It can be shown that the 
normalized eigenvector v = (sp ,yo) (|v| — 1) associated with Ao points in the direction of the largest amount of 
polarization [6l| . One can define an elliptical component of polarization Pg as 



ft = ^, (D5) 



14 In practice we do not use the gravitational wave amplitudes, but the real and imaginary components of <E I 4, in order to avoid problems 
with the integration constants and to sidestep the memory effect discussed in Section III AI This should have no influence on the final 
results. 
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where 

X = max ^Rc[e la x ] 2 + Re[e la y } 2 . (D6) 

a 

The quantity Pe — 1 for circular polarization, and Pg = for linear polarization. For illustration, consider the 
waveform h+ = sin i,/i x = cost. This is a good approximation to a typical inspiral waveform at large orbital 
separation, as viewed from the normal to the orbital plane, and it is obviously circularly polarized. For this waveform 
we have 

c =0?)' (° 7 ) 

and v = (— i/V2, l/y/2). This implies X = l/v2 and Pe = 1, as expected. For h x = (linear polarization) we would 
have X = 1 and Pe — 0. Thus Pe is a good indicator of the degree of circular or linear polarization. We can also 
define a polarization strength as follows: 

P S = 1-^. (D8) 
Ao 

Ps = 1 means that the waveform is entirely linearly polarized or circularly polarized (there is only one polarization 
component), and Ps = means that the two polarization states have comparable magnitude. 




FIG. 28: Top: |Re(Afr tp22)\ and |Im(Mr ^/>22)| for q — 2.0, D8 (top). Bottom: the degree of elliptic polarization Pe for the 
1 = 2 waveform as viewed from the normal to the orbital plane. 



In Fig. [28]we show the result of computing Pe using the dominant (/ = \m\ — 2) component of a binary black hole 
merger waveform with q = 2.0. This plot clearly shows that the polarization is circular for both inspiral and ringdown, 
with the exception of the unphysical portions of the wave: the initial data burst and the final, noise-dominated part 
of the ringdown waveform. 
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